SRM_Workshop_2022.ipynb

Kernel: Python 3

Interactive mapping and analysis of geospatial big data using geemap and Google Earth Engine

This notebook was developed for the geemap workshop at the Society for Range Management (SRM) 2022 Annual Meeting.

Authors: Qiusheng Wu

Link to this notebook: https://gishub.org/SRM

Introduction

Description

Google Earth Engine (GEE) is a cloud computing platform with a multi-petabyte catalog of satellite imagery and geospatial datasets. It enables scientists, researchers, and developers to analyze and visualize changes on the Earth’s surface. The geemap Python package provides GEE users with an intuitive interface to manipulate, analyze, and visualize geospatial big data interactively in a Jupyter-based environment. The topics to be covered in this workshop include:

1. Introducing geemap

2. Creating interactive maps

3. Searching GEE data catalog

4. Data Visualization

5. Data Analysis

6. Data Export

This workshop is intended for scientific programmers, data scientists, geospatial analysts, and concerned citizens of Earth. The attendees are expected to have a basic understanding of Python and the Jupyter ecosystem. Familiarity with Earth science and geospatial datasets is useful but not required.

Useful links

• Streamlit Geospatial

• Google Earth Engine

• geemap.org

• Google Earth Engine and geemap Python Tutorials (55 videos with a total length of 15 hours)

• Spatial Data Management with Google Earth Engine (19 videos with a total length of 9 hours)

• Ask geemap questions on GitHub

Prerequisite

• A Google Earth Engine account. Sigh up here if needed.

• Miniconda or Anaconda has been installed on your computer.

Set up a conda environment

conda create -n gee python=3.9
conda activate gee
conda install geopandas
conda install mamba -c conda-forge
mamba install geemap localtileserver -c conda-forge

Launch JupyterLab

Open Anaconda Prompt or the Terminal and enter the following commands.

conda activate gee
jupyter lab

Geemap basics

Install geemap

import subprocess

try:
    import geemap
except ImportError:
    print("Installing geemap ...")
    subprocess.check_call(["python", "-m", "pip", "install", "geemap"])

Import libraries

import os
import ee
import geemap

Create an interactive map

Map = geemap.Map()
Map

Customize the default map

You can specify the center(lat, lon) and zoom for the default map. You can also set the visibility of the controls.

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

Add basemaps

Map = geemap.Map()
Map.add_basemap("HYBRID")
Map.add_basemap("OpenTopoMap")
Map

Change basemaps interactively

Map = geemap.Map()
Map

Use drawing tools

Map = geemap.Map()
Map

# Map.user_roi.getInfo()

# Map.user_rois.getInfo()

Convert GEE JavaScript to Python

https://developers.google.com/earth-engine/guides/image_visualization

js_snippet = """
// Load an image.
var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318');

// Define the visualization parameters.
var 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');

"""

geemap.js_snippet_to_py(
    js_snippet, add_new_cell=True, import_ee=True, import_geemap=True, show_map=True
)

import ee
import geemap

Map = geemap.Map()

# Load an image.
image = ee.Image("LANDSAT/LC08/C01/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")
Map

You can also convert GEE JavaScript to Python without coding.

Map = geemap.Map()
Map

Earth Engine datasets

Load Earth Engine datasets

More GEE datasets can be found at https://developers.google.com/earth-engine/datasets/

Map = geemap.Map()

# Add Earth Engine datasets
dem = ee.Image("USGS/SRTMGL1_003")
landsat7 = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")
states = ee.FeatureCollection("TIGER/2018/States")

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

# Add Earth Engine layers to Map
Map.addLayer(dem, vis_params, "SRTM DEM", True, 0.5)
Map.addLayer(
    landsat7,
    {"bands": ["B4", "B3", "B2"], "min": 20, "max": 200, "gamma": 1.5},
    "Landsat 7",
)
Map.addLayer(states, {}, "US States")

Map

Search the Earth Engine Data Catalog

Map = geemap.Map()
Map

dem = ee.Image("CGIAR/SRTM90_V4")
Map.addLayer(dem, {}, "CGIAR/SRTM90_V4")

hillshade = ee.Terrain.hillshade(dem)
Map.addLayer(hillshade, {}, "Hillshade")

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

Map.addLayer(dem, vis_params, "DEM", True, 0.5)

Use the datasets module

from geemap.datasets import DATA

Map = geemap.Map()

dem = ee.Image(DATA.USGS_SRTMGL1_003)

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

Map.addLayer(dem, vis_params, "SRTM DEM")
Map

Use the Inspector tool

Map = geemap.Map()

# Add Earth Engine datasets
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")

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

# Add Earth Engine layers to Map
Map.addLayer(dem, vis_params, "SRTM DEM", True, 0.5)
Map.addLayer(
    landsat7,
    {"bands": ["B4", "B3", "B2"], "min": 20, "max": 180, "gamma": 1.5},
    "Landsat 7",
)
Map.addLayer(states, {}, "US States")

Map

Data visualization

Use the Plotting tool

Map = geemap.Map()

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

Change layer opacity

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

dem = ee.Image("USGS/SRTMGL1_003")
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", True, 1)
Map.addLayer(states, {}, "US States", True)

Map

Visualize raster data

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

# Add Earth Engine dataset
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"],
}

Map.addLayer(dem, vis_params, "SRTM DEM", True, 1)
Map.addLayer(
    landsat7,
    {"bands": ["B4", "B3", "B2"], "min": 20, "max": 200, "gamma": 2},
    "Landsat 7",
)
Map

Visualize vector data

Map = geemap.Map()

states = ee.FeatureCollection("TIGER/2018/States")
Map.addLayer(states, {}, "US States")
Map

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
)

Add a legend

legends = geemap.builtin_legends
for legend in legends:
    print(legend)

Map = geemap.Map()
Map.add_basemap("HYBRID")

dataset = ee.ImageCollection("USGS/NLCD_RELEASES/2016_REL")
nlcd2016 = dataset.filter(ee.Filter.eq("system:index", "2016")).first()
landcover = nlcd2016.select("landcover")

Map.addLayer(landcover, {}, "NLCD Land Cover 2016")
Map.add_legend(builtin_legend="NLCD")
Map

Map = geemap.Map()

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",
}

dataset = ee.ImageCollection("USGS/NLCD_RELEASES/2016_REL")
nlcd2016 = dataset.filter(ee.Filter.eq("system:index", "2016")).first()
landcover = nlcd2016.select("landcover")

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

Add a colorbar

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

colors = vis_params["palette"]
vmin = vis_params["min"]
vmax = vis_params["max"]

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

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

Map.add_colorbar(
    vis_params,
    label="Elevation (m)",
    layer_name="SRTM DEM",
    orientation="vertical",
    transparent_bg=True,
)

Map.add_colorbar(
    vis_params,
    label="Elevation (m)",
    layer_name="SRTM DEM",
    orientation="vertical",
    transparent_bg=True,
    discrete=True,
)

Create a split-panel map

Map = geemap.Map()
Map.split_map(left_layer="HYBRID", right_layer="TERRAIN")
Map

Map = geemap.Map(center=(40, -100), zoom=4)
Map.split_map(
    left_layer="NLCD 2016 CONUS Land Cover", right_layer="NLCD 2001 CONUS Land Cover"
)
Map

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

dataset = ee.ImageCollection("USGS/NLCD_RELEASES/2016_REL")
nlcd_2001 = (
    dataset.filter(ee.Filter.eq("system:index", "2001")).first().select("landcover")
)
nlcd_2016 = (
    dataset.filter(ee.Filter.eq("system:index", "2016")).first().select("landcover")
)

left_layer = geemap.ee_tile_layer(nlcd_2001, {}, "NLCD 2001")
right_layer = geemap.ee_tile_layer(nlcd_2016, {}, "NLCD 2016")

Map.split_map(left_layer, right_layer)
Map

Create linked maps

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="400px",
    center=[38.4151, 21.2712],
    zoom=12,
    ee_objects=[image],
    vis_params=vis_params,
    labels=labels,
    label_position="topright",
)

Use timeseries inspector

Using data from the Rangeland Analysis Platform.

• Vegetation cover dataset: ee.ImageCollection('projects/rangeland-analysis-platform/vegetation-cover-v3')

• Rangeland production dataset: ee.ImageCollection('projects/rangeland-analysis-platform/npp-partitioned-v3')

import geemap.colormaps as cm

Map = geemap.Map()

collection = ee.ImageCollection(
    "projects/rangeland-analysis-platform/vegetation-cover-v3"
).select("TRE")

years = collection.aggregate_array("year").getInfo()
years

names = [str(i) for i in years]
names

palette = cm.palettes.ndvi
palette

vis_params = {"min": 0, "max": 50, "palette": palette}

# Map.addLayer(collection.first(), {}, "First image")
Map.ts_inspector(
    left_ts=collection,
    right_ts=collection,
    left_names=names,
    right_names=names,
    left_vis=vis_params,
    right_vis=vis_params,
)
Map

Data analysis

Descriptive statistics

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

image.propertyNames().getInfo()

image.get("CLOUD_COVER").getInfo()

props = geemap.image_props(image)
props.getInfo()

stats = geemap.image_stats(image, scale=90)
stats.getInfo()

Zonal statistics

Map = geemap.Map()

# Add Earth Engine dataset
dem = ee.Image("USGS/SRTMGL1_003")

# Set visualization parameters.
dem_vis = {
    "min": 0,
    "max": 4000,
    "palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}

# Add Earth Engine DEM to map
Map.addLayer(dem, dem_vis, "SRTM DEM")

# Add Landsat data to map
landsat = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")

landsat_vis = {"bands": ["B4", "B3", "B2"], "gamma": 1.4}
Map.addLayer(landsat, landsat_vis, "LE7_TOA_5YEAR/1999_2003")

states = ee.FeatureCollection("TIGER/2018/States")
Map.addLayer(states, {}, "US States")
Map

out_dem_stats = "dem_stats.csv"

# Allowed output formats: csv, shp, json, kml, kmz
# Allowed statistics type: MEAN, MAXIMUM, MINIMUM, MEDIAN, STD, MIN_MAX, VARIANCE, SUM
geemap.zonal_stats(dem, states, out_dem_stats, stat_type="MEAN", scale=1000)

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

Zonal statistics by group

Map = geemap.Map()

dataset = ee.ImageCollection("USGS/NLCD_RELEASES/2016_REL")
nlcd2016 = dataset.filter(ee.Filter.eq("system:index", "2016")).first()
landcover = nlcd2016.select("landcover")
Map.addLayer(landcover, {}, "NLCD 2016")

states = ee.FeatureCollection("TIGER/2018/States")
Map.addLayer(states, {}, "US States")
Map.add_legend(builtin_legend="NLCD")
Map

nlcd_stats = "nlcd_stats.csv"

# statistics_type can be either 'SUM' or 'PERCENTAGE'
# denominator can be used to convert square meters to other areal units, such as square kilimeters
geemap.zonal_stats_by_group(
    landcover,
    states,
    nlcd_stats,
    stat_type="SUM",
    denominator=1000000,
    decimal_places=2,
)

WhiteboxTools

https://www.whiteboxgeo.com/manual/wbt_book/intro.html

import whiteboxgui

whiteboxgui.show()

whiteboxgui.show(tree=True)

Map = geemap.Map()
Map

Data export

Export ee.Image

Map = geemap.Map()

image = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")

landsat_vis = {"bands": ["B4", "B3", "B2"], "gamma": 1.4}
Map.addLayer(image, landsat_vis, "LE7_TOA_5YEAR/1999_2003", True, 1)

Map

# Draw any shapes on the map using the Drawing tools before executing this code block
roi = Map.user_roi

if roi is None:
    roi = ee.Geometry.Polygon(
        [
            [
                [-115.413031, 35.889467],
                [-115.413031, 36.543157],
                [-114.034328, 36.543157],
                [-114.034328, 35.889467],
                [-115.413031, 35.889467],
            ]
        ]
    )

filename = "landsat.tif"

Exporting all bands as one single image

image = image.clip(roi).unmask()
geemap.ee_export_image(
    image, filename=filename, scale=90, region=roi, file_per_band=False
)

Exporting each band as one image

geemap.ee_export_image(
    image, filename=filename, scale=90, region=roi, file_per_band=True
)

Export an image to Google Drive¶

# geemap.ee_export_image_to_drive(image, description='landsat', folder='export', region=roi, scale=30)

Export ee.ImageCollection

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

collection.aggregate_array("system:index").getInfo()

out_dir = os.getcwd()
geemap.ee_export_image_collection(collection, out_dir=out_dir)

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

Extract pixels as a numpy array

import matplotlib.pyplot as plt

img = ee.Image("LANDSAT/LC08/C01/T1_SR/LC08_038029_20180810").select(["B4", "B5", "B6"])

aoi = ee.Geometry.Polygon(
    [[[-110.8, 44.7], [-110.8, 44.6], [-110.6, 44.6], [-110.6, 44.7]]], None, False
)

rgb_img = geemap.ee_to_numpy(img, region=aoi)
print(rgb_img.shape)

rgb_img_test = (255 * ((rgb_img[:, :, 0:3] - 100) / 3500)).astype("uint8")
plt.imshow(rgb_img_test)
plt.show()

Export pixel values to points

Map = geemap.Map()

# Add Earth Engine dataset
dem = ee.Image("USGS/SRTMGL1_003")
landsat7 = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")

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

# Add Earth Engine layers to Map
Map.addLayer(
    landsat7, {"bands": ["B4", "B3", "B2"], "min": 20, "max": 200}, "Landsat 7"
)
Map.addLayer(dem, vis_params, "SRTM DEM", True, 1)
Map

Download sample data

work_dir = os.getcwd()
in_shp = os.path.join(work_dir, "us_cities.shp")
if not os.path.exists(in_shp):
    data_url = "https://github.com/giswqs/data/raw/main/us/us_cities.zip"
    geemap.download_from_url(data_url, out_dir=work_dir)

in_fc = geemap.shp_to_ee(in_shp)
Map.addLayer(in_fc, {}, "Cities")

Export pixel values as a shapefile

out_shp = os.path.join(work_dir, "dem.shp")
geemap.extract_values_to_points(in_fc, dem, out_shp)

Export pixel values as a csv

out_csv = os.path.join(work_dir, "landsat.csv")
geemap.extract_values_to_points(in_fc, landsat7, out_csv)

Export ee.FeatureCollection

Map = geemap.Map()

fc = ee.FeatureCollection("users/giswqs/public/countries")
Map.addLayer(fc, {}, "Countries")
Map

out_shp = "countries.shp"

geemap.ee_to_shp(fc, filename=out_shp)

out_csv = os.path.join(out_dir, "countries.csv")
geemap.ee_export_vector(fc, filename=out_csv)

out_kml = os.path.join(out_dir, "countries.kml")
geemap.ee_export_vector(fc, filename=out_kml)

# geemap.ee_export_vector_to_drive(fc, description="countries", folder="export", file_format="shp")