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.

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:

Introducing geemap
Creating interactive maps
Searching GEE data catalog
Data Visualization
Data Analysis
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¶

In [ ]:

import subprocess

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

Import libraries¶

In [ ]:

import os
import ee
import geemap

Create an interactive map¶

In [ ]:

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.

In [ ]:

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

Add basemaps¶

In [ ]:

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

Change basemaps interactively¶

In [ ]:

Map = geemap.Map()
Map

Use drawing tools¶

In [ ]:

Map = geemap.Map()
Map

In [ ]:

# Map.user_roi.getInfo()

In [ ]:

# Map.user_rois.getInfo()

Convert GEE JavaScript to Python¶

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

In [ ]:

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');

"""

In [ ]:

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

In [ ]:

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.

In [ ]:

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/

In [ ]:

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¶

In [ ]:

Map = geemap.Map()
Map

In [ ]:

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

In [ ]:

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

In [ ]:

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¶

In [ ]:

from geemap.datasets import DATA

In [ ]:

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¶

In [ ]:

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¶

In [ ]:

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¶

In [ ]:

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¶

In [ ]:

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¶

In [ ]:

Map = geemap.Map()

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

In [ ]:

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¶

In [ ]:

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

In [ ]:

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

In [ ]:

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¶

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

Create a split-panel map¶

In [ ]:

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

In [ ]:

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

In [ ]:

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¶

In [ ]:

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

In [ ]:

import geemap.colormaps as cm

In [ ]:

Map = geemap.Map()

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

palette = cm.palettes.ndvi
palette

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

image.propertyNames().getInfo()

In [ ]:

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

In [ ]:

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

In [ ]:

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

Zonal statistics¶

In [ ]:

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

In [ ]:

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)

In [ ]:

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

Zonal statistics by group¶

In [ ]:

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

In [ ]:

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

In [ ]:

import whiteboxgui

In [ ]:

whiteboxgui.show()

In [ ]:

whiteboxgui.show(tree=True)

In [ ]:

Map = geemap.Map()
Map

Data export¶

Export ee.Image¶

In [ ]:

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

In [ ]:

# 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],
            ]
        ]
    )

In [ ]:

filename = "landsat.tif"

Exporting all bands as one single image

In [ ]:

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

In [ ]:

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

Export an image to Google Drive¶

In [ ]:

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

Export ee.ImageCollection¶

In [ ]:

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

In [ ]:

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

In [ ]:

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

In [ ]:

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

Extract pixels as a numpy array¶

In [ ]:

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)

In [ ]:

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

Export pixel values to points¶

In [ ]:

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

In [ ]:

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 [ ]:

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

Export pixel values as a shapefile

In [ ]:

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

In [ ]:

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

Export ee.FeatureCollection¶

In [ ]:

Map = geemap.Map()

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

In [ ]:

out_shp = "countries.shp"

In [ ]:

geemap.ee_to_shp(fc, filename=out_shp)

In [ ]:

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

In [ ]:

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

In [ ]:

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