Geospatial Cloud Computing with the GEE Python API - Part 2

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

Introduction¶

This notebook contains the materials for the first 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 second part will cover the following topics:

Processing of vector data (shapefiles, json, conversion from one format to another)
Processing of raster data: extract pixel value, raster calculator, zonal statistics etc.
Working with local geospatial data in Geemap
Accessing Cloud Optimized GeoTIFF
Exporting EE Image and Feature data
Creating timelapse animations using Landsat or Sentinel 2 for Alaska
Time series analysis: Forest cover change for a test site in Alaska (e.g. Bonanza Creek LTER or Caribou-Poker Creeks Research Watershed)

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

In [ ]:

# %pip install geemap pygis

Import libraries¶

In [ ]:

import ee
import geemap

In [ ]:

geemap.ee_initialize()

Processing of vector data¶

From GeoJSON¶

In [ ]:

in_geojson = "https://github.com/gee-community/geemap/blob/master/examples/data/countries.geojson"
m = geemap.Map()
fc = geemap.geojson_to_ee(in_geojson)
m.add_layer(fc.style(**{"color": "ff0000", "fillColor": "00000000"}), {}, "Countries")
m

From Shapefile¶

In [ ]:

url = "https://github.com/gee-community/geemap/blob/master/examples/data/countries.zip"
geemap.download_file(url, overwrite=True)

In [ ]:

in_shp = "countries.shp"
fc = geemap.shp_to_ee(in_shp)

In [ ]:

m = geemap.Map()
m.add_layer(fc, {}, "Countries")
m

From GeoDataFrame¶

In [ ]:

import geopandas as gpd

gdf = gpd.read_file(in_shp)
fc = geemap.gdf_to_ee(gdf)

In [ ]:

m = geemap.Map()
m.add_layer(fc, {}, "Countries")
m

To GeoJSON¶

In [ ]:

m = geemap.Map()
states = ee.FeatureCollection("TIGER/2018/States")
fc = states.filter(ee.Filter.eq("NAME", "Alaska"))
m.add_layer(fc, {}, "Alaska")
m.center_object(fc, 4)
m

In [ ]:

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

To Shapefile¶

In [ ]:

geemap.ee_to_shp(fc, filename="Alaska.shp")

To GeoDataFrame¶

In [ ]:

gdf = geemap.ee_to_gdf(fc)
gdf

In [ ]:

gdf.explore()

To DataFrame¶

In [ ]:

df = geemap.ee_to_df(fc)
df

To CSV¶

In [ ]:

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

Processing of raster data¶

Extract pixel values¶

Extracting values to points¶

In [ ]:

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

dem = ee.Image("USGS/SRTMGL1_003")
landsat7 = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")

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 7",
)
m.add_layer(dem, vis_params, "SRTM DEM", True, 1)
m

In [ ]:

in_shp = "us_cities.shp"
url = "https://github.com/giswqs/data/raw/main/us/us_cities.zip"
geemap.download_file(url)

In [ ]:

in_fc = geemap.shp_to_ee(in_shp)
m.add_layer(in_fc, {}, "Cities")

In [ ]:

geemap.extract_values_to_points(in_fc, dem, out_fc="dem.shp")

In [ ]:

geemap.shp_to_gdf("dem.shp")

In [ ]:

geemap.extract_values_to_points(in_fc, landsat7, "landsat.csv")

In [ ]:

geemap.csv_to_df("landsat.csv")

Extracting pixel values along a transect¶

In [ ]:

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

image = ee.Image("USGS/SRTMGL1_003")
vis_params = {
    "min": 0,
    "max": 4000,
    "palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}
m.add_layer(image, vis_params, "SRTM DEM", True, 0.5)
m

In [ ]:

line = m.user_roi
if line is None:
    line = ee.Geometry.LineString(
        [[-120.2232, 36.3148], [-118.9269, 36.7121], [-117.2022, 36.7562]]
    )
    m.add_layer(line, {}, "ROI")
m.centerObject(line)

In [ ]:

reducer = "mean"
transect = geemap.extract_transect(
    image, line, n_segments=100, reducer=reducer, to_pandas=True
)
transect

In [ ]:

geemap.line_chart(
    data=transect,
    x="distance",
    y="mean",
    markers=True,
    x_label="Distance (m)",
    y_label="Elevation (m)",
    height=400,
)

In [ ]:

transect.to_csv("transect.csv")

Zonal statistics¶

Zonal statistics with an image and a feature collection¶

In [ ]:

m = 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"],
}
m.add_layer(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}
m.add_layer(landsat, landsat_vis, "Landsat", False)

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

In [ ]:

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

In [ ]:

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

Zonal statistics by group¶

In [ ]:

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

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

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

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

In [ ]:

nlcd_stats = "nlcd_stats.csv"

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

In [ ]:

nlcd_stats = "nlcd_stats_pct.csv"

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

Zonal statistics with two images¶

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

Map algebra¶

In [ ]:

m = geemap.Map()

# Load a 5-year Landsat 7 composite 1999-2003.
landsat_1999 = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")

# Compute NDVI.
ndvi_1999 = (
    landsat_1999.select("B4")
    .subtract(landsat_1999.select("B3"))
    .divide(landsat_1999.select("B4").add(landsat_1999.select("B3")))
)

vis = {"min": 0, "max": 1, "palette": "ndvi"}
m.add_layer(ndvi_1999, vis, "NDVI")
m.add_colorbar(vis, label="NDVI")
m

In [ ]:

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

# Compute the EVI using an expression.
evi = image.expression(
    "2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))",
    {
        "NIR": image.select("B5"),
        "RED": image.select("B4"),
        "BLUE": image.select("B2"),
    },
)

# Define a map centered on San Francisco Bay.
m = geemap.Map(center=[37.4675, -122.1363], zoom=9)

vis = {"min": 0, "max": 1, "palette": "ndvi"}
m.add_layer(evi, vis, "EVI")
m.add_colorbar(vis, label="EVI")
m

Exercise 1 - Zonal statistics¶

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

Working with local geospatial data¶

Raster data¶

Single-band raster¶

In [ ]:

url = "https://github.com/giswqs/data/raw/main/raster/srtm90.tif"
filename = "dem.tif"
geemap.download_file(url, filename)

Multi-band raster¶

In [ ]:

m = geemap.Map()
m.add_raster(filename, cmap="terrain", layer_name="DEM")
vis_params = {"min": 0, "max": 4000, "palette": "terrain"}
m.add_colorbar(vis_params, label="Elevation (m)")
m

In [ ]:

url = "https://github.com/giswqs/data/raw/main/raster/cog.tif"
filename = "cog.tif"
geemap.download_file(url, filename)

In [ ]:

m = geemap.Map()
m.add_raster(filename, indexes=[4, 1, 2], layer_name="False color")
m

Vector data¶

GeoJSON¶

In [ ]:

in_geojson = (
    "https://github.com/opengeos/datasets/releases/download/vector/cables.geojson"
)
m = geemap.Map()
m.add_geojson(in_geojson, layer_name="Cable lines", info_mode="on_hover")
m

In [ ]:

m = geemap.Map()
m.add_basemap("CartoDB.DarkMatter")
callback = lambda feat: {"color": feat["properties"]["color"], "weight": 2}
m.add_geojson(in_geojson, layer_name="Cable lines", style_callback=callback)
m

In [ ]:

url = "https://github.com/opengeos/datasets/releases/download/world/countries.geojson"
m = geemap.Map()
m.add_geojson(
    url, layer_name="Countries", fill_colors=["red", "yellow", "green", "orange"]
)
m

In [ ]:

import random

m = geemap.Map()


def random_color(feature):
    return {
        "color": "black",
        "weight": 3,
        "fillColor": random.choice(["red", "yellow", "green", "orange"]),
    }


m.add_geojson(url, layer_name="Countries", style_callback=random_color)
m

In [ ]:

m = geemap.Map()

style = {
    "stroke": True,
    "color": "#0000ff",
    "weight": 2,
    "opacity": 1,
    "fill": True,
    "fillColor": "#0000ff",
    "fillOpacity": 0.1,
}

hover_style = {"fillOpacity": 0.7}

m.add_geojson(url, layer_name="Countries", style=style, hover_style=hover_style)
m

Shapefile¶

In [ ]:

url = "https://github.com/opengeos/datasets/releases/download/world/countries.zip"
geemap.download_file(url, overwrite=True)

In [ ]:

m = geemap.Map()
in_shp = "countries.shp"
m.add_shp(in_shp, layer_name="Countries")
m

GeoDataFrame¶

In [ ]:

import geopandas as gpd

m = geemap.Map(center=[40, -100], zoom=4)
gdf = gpd.read_file("countries.shp")
m.add_gdf(gdf, layer_name="Countries")
m

GeoPackage¶

In [ ]:

m = geemap.Map()
data = "https://github.com/opengeos/datasets/releases/download/world/countries.gpkg"
m.add_vector(data, layer_name="Countries")
m

CSV to vector¶

In [ ]:

data = "https://github.com/gee-community/geemap/blob/master/examples/data/us_cities.csv"
geemap.csv_to_df(data)

In [ ]:

geemap.csv_to_geojson(
    data, "cities.geojson", latitude="latitude", longitude="longitude"
)

In [ ]:

geemap.csv_to_shp(data, "cities.shp", latitude="latitude", longitude="longitude")

In [ ]:

geemap.csv_to_vector(data, "cities.gpkg", latitude="latitude", longitude="longitude")

In [ ]:

gdf = geemap.csv_to_gdf(data, latitude="latitude", longitude="longitude")
gdf

In [ ]:

cities = (
    "https://github.com/gee-community/geemap/blob/master/examples/data/us_cities.csv"
)
m = geemap.Map(center=[40, -100], zoom=4)
m.add_points_from_xy(cities, x="longitude", y="latitude")
m

In [ ]:

regions = "https://github.com/gee-community/geemap/blob/master/examples/data/us_regions.geojson"

In [ ]:

m = geemap.Map(center=[40, -100], zoom=4)
m.add_geojson(regions, layer_name="US Regions")
m.add_points_from_xy(
    cities,
    x="longitude",
    y="latitude",
    layer_name="US Cities",
    color_column="region",
    icon_names=["gear", "map", "leaf", "globe"],
    spin=True,
    add_legend=True,
)
m

In [ ]:

m = geemap.Map(center=[40, -100], zoom=4)
m.add_circle_markers_from_xy(
    data,
    x="longitude",
    y="latitude",
    radius=8,
    color="blue",
    fill_color="black",
    fill_opacity=0.5,
)
m

Accessing Cloud Optimized GeoTIFFs¶

COG¶

In [ ]:

m = geemap.Map()
url = "https://opendata.digitalglobe.com/events/california-fire-2020/pre-event/2018-02-16/pine-gulch-fire20/1030010076004E00.tif"
m.add_cog_layer(url, name="Fire (pre-event)")
m

In [ ]:

geemap.cog_center(url)

In [ ]:

geemap.cog_bands(url)

In [ ]:

url2 = "https://opendata.digitalglobe.com/events/california-fire-2020/post-event/2020-08-14/pine-gulch-fire20/10300100AAC8DD00.tif"
m.add_cog_layer(url2, name="Fire (post-event)")

In [ ]:

m = geemap.Map(center=[39.4568, -108.5107], zoom=12)
m.split_map(left_layer=url2, right_layer=url)
m

SpatioTemporal Asset Catalog (STAC)¶

In [ ]:

url = "https://tinyurl.com/22vptbws"

In [ ]:

geemap.stac_bounds(url)

In [ ]:

geemap.stac_center(url)

In [ ]:

geemap.stac_bands(url)

In [ ]:

m = geemap.Map()
m.add_stac_layer(url, bands=["pan"], name="Panchromatic")
m.add_stac_layer(url, bands=["B3", "B2", "B1"], name="False color")
m

Exporting Earth Engine data¶

Exporting images¶

Add a Landsat image to the map.

In [ ]:

m = 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]}

m.center_object(image)
m.add_layer(image, vis_params, "Landsat")
m

Add a rectangle to the map.

In [ ]:

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

To local drive

In [ ]:

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

Check image projection.

In [ ]:

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

In [ ]:

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

Specify region, crs, and crs_transform.

In [ ]:

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

To Google Drive

In [ ]:

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

In [ ]:

geemap.download_ee_image(image, "landsat.tif", scale=90)

Exporting image collections¶

In [ ]:

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

In [ ]:

collection.aggregate_array("system:index")

To local drive

In [ ]:

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

To Google Drive

In [ ]:

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

Exporting feature collections¶

In [ ]:

m = geemap.Map()
states = ee.FeatureCollection("TIGER/2018/States")
fc = states.filter(ee.Filter.eq("NAME", "Alaska"))
m.add_layer(fc, {}, "Alaska")
m.center_object(fc, 4)
m

To local drive

In [ ]:

geemap.ee_to_shp(fc, filename="Alaska.shp")

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

gdf = geemap.ee_to_gdf(fc)
gdf

In [ ]:

df = geemap.ee_to_df(fc)
df

To Google Drive

In [ ]:

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

Creating timelapse animations¶

Landsat timelapse¶

In [ ]:

m = geemap.Map(center=[64.838721, -147.763366], zoom=11)
m

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.

In [ ]:

roi = m.user_roi
if roi is None:
    roi = ee.Geometry.BBox(-147.9701, 64.7733, -147.5849, 64.8717)
    m.add_layer(roi)
    m.center_object(roi)

In [ ]:

timelapse = geemap.landsat_timelapse(
    roi,
    out_gif="Fairbanks.gif",
    start_year=2000,
    end_year=2023,
    start_date="06-01",
    end_date="09-01",
    bands=["SWIR1", "NIR", "Red"],
    frames_per_second=5,
    title="Landsat Timelapse",
    progress_bar_color="blue",
    mp4=True,
)
geemap.show_image(timelapse)

In [ ]:

m = geemap.Map(center=[64.838721, -147.763366], zoom=11)
m.add_gui("timelapse")
m

In [ ]:

m = geemap.Map()
roi = ee.Geometry.BBox(-115.5541, 35.8044, -113.9035, 36.5581)
m.add_layer(roi)
m.center_object(roi)
m

In [ ]:

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

In [ ]:

m = geemap.Map()
roi = ee.Geometry.BBox(113.8252, 22.1988, 114.0851, 22.3497)
m.add_layer(roi)
m.center_object(roi)
m

In [ ]:

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¶

In [ ]:

m = geemap.Map(center=[64.838721, -147.763366], zoom=11)
m

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.

In [ ]:

roi = m.user_roi
if roi is None:
    roi = ee.Geometry.BBox(-147.9701, 64.7733, -147.5849, 64.8717)
    m.add_layer(roi)
    m.center_object(roi)

In [ ]:

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

MODIS¶

MODIS vegetation indices

In [ ]:

m = geemap.Map()
m

In [ ]:

roi = m.user_roi
if roi is None:
    roi = ee.Geometry.BBox(-18.6983, -36.1630, 52.2293, 38.1446)
    m.add_layer(roi)
    m.center_object(roi)

In [ ]:

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

In [ ]:

m = geemap.Map()
m

In [ ]:

roi = m.user_roi
if roi is None:
    roi = ee.Geometry.BBox(-171.21, -57.13, 177.53, 79.99)
    m.add_layer(roi)
    m.center_object(roi)

In [ ]:

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¶

In [ ]:

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"

In [ ]:

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

In [ ]:

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"

In [ ]:

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

In [ ]:

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"

In [ ]:

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

Exercise 2 - 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.

In [ ]:

Time series analysis¶

Visualizing forest cover¶

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

In [ ]:

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

Select the imagery for 2000.

In [ ]:

m = geemap.Map()
first_bands = ["first_b50", "first_b40", "first_b30"]
first_image = dataset.select(first_bands)
m.add_layer(first_image, {"bands": first_bands, "gamma": 1.5}, "Landsat 2000")

Select the imagery for 2022.

In [ ]:

last_bands = ["last_b50", "last_b40", "last_b30"]
last_image = dataset.select(last_bands)
m.add_layer(last_image, {"bands": last_bands, "gamma": 1.5}, "Landsat 2022")

Select the tree cover imagery for 2000.

In [ ]:

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

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

In [ ]:

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

Visualizing forest gain and loss¶

Visualize forest loss.

In [ ]:

m = geemap.Map(center=[64.864983, -147.840441], zoom=4)
m.add_basemap("Esri.WorldImagery")
treeloss_year = dataset.select(["lossyear"])
treeLossVisParam = {"min": 0, "max": 22, "palette": ["yellow", "red"]}
layer_name = "Tree loss year"
m.add_layer(treeloss_year, treeLossVisParam, layer_name)
m.add_colorbar(treeLossVisParam, label=layer_name, layer_name=layer_name)
m.add("layer_manager")
m

Compare forest loss and gain.

In [ ]:

m = geemap.Map(center=[64.864983, -147.840441], zoom=4)
m.add_basemap("Esri.WorldImagery")
treeloss = dataset.select(["loss"]).selfMask()
treegain = dataset.select(["gain"]).selfMask()
m.add_layer(treeloss, {"palette": "red"}, "Tree loss")
m.add_layer(treegain, {"palette": "yellow"}, "Tree gain")
m.add("layer_manager")
m

Calculating forest cover change¶

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

Add a county boundary layer to the map.

In [ ]:

counties = ee.FeatureCollection("TIGER/2018/Counties").filter(
    ee.Filter.eq("STATEFP", "02")
)
df = geemap.ee_to_df(counties)
df.head()

In [ ]:

style = {"color": "0000FFFF", "fillColor": "00000000"}
m.add_layer(counties, {}, "Counties Vector", False)
m.add_layer(counties.style(**style), {}, "Counties Raster")
m

Compute zonal statistics by county.

In [ ]:

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

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

In [ ]:

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

In [ ]:

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

Calculate the forest loss area by county.

In [ ]:

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

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

In [ ]:

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

In [ ]:

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

Exercise 3 - 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:

In [ ]: