Alaska_2024_Part2.ipynb

Kernel: Python 3

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

# %pip install geemap pygis

Import libraries

import ee
import geemap

geemap.ee_initialize()

Processing of vector data

From GeoJSON

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

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

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

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

From GeoDataFrame

import geopandas as gpd

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

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

To GeoJSON

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

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

To Shapefile

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

To GeoDataFrame

gdf = geemap.ee_to_gdf(fc)
gdf

gdf.explore()

To DataFrame

df = geemap.ee_to_df(fc)
df

To CSV

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

Processing of raster data

Extract pixel values

Extracting values to points

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_shp = "us_cities.shp"
url = "https://github.com/giswqs/data/raw/main/us/us_cities.zip"
geemap.download_file(url)

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

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

geemap.shp_to_gdf("dem.shp")

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

geemap.csv_to_df("landsat.csv")

Extracting pixel values along a transect

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

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)

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

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

transect.to_csv("transect.csv")

Zonal statistics

Zonal statistics with an image and a feature collection

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

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

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

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

nlcd_stats = "nlcd_stats.csv"

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

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

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

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

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

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

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

Map algebra

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

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

• ee.FeatureCollection("TIGER/2018/States")

• ee.ImageCollection("NOAA/GFS0P25")

Working with local geospatial data

Raster data

Single-band raster

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

Multi-band raster

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

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

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

Vector data

GeoJSON

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

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

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

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

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

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

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

GeoDataFrame

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

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

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

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

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

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

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

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

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

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

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

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

geemap.cog_center(url)

geemap.cog_bands(url)

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

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

SpatioTemporal Asset Catalog (STAC)

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

geemap.stac_bounds(url)

geemap.stac_center(url)

geemap.stac_bands(url)

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.

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.

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

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
)

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

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

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

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

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

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

geemap.ee_to_csv(fc, filename="Alaska.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="Alaska", fileFormat="SHP", folder="export"
)

Creating timelapse animations

Landsat timelapse

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.

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)

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)

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

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

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)

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

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

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.

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)

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

m = geemap.Map()
m

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)

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

m = geemap.Map()
m

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)

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

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

Time series analysis

Visualizing forest cover

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.

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.

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.

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

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.

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.

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.

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

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.

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.

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.

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.

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.

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.

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:

• ee.FeatureCollection("TIGER/2018/States")

• ee.Image("UMD/hansen/global_forest_change_2022_v1_10")