Geospatial Cloud Computing with the GEE Python API - Part 1
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.
The workshop is divided into three parts. The first part will cover the following topics:
Google Earth Engine is a cloud-computing platform for scientific analysis and visualization of geospatial datasets. It 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.
Currently, ipywidgets does not work well with Colab dark theme. Some of the geemap widgets may not display properly in Colab dark theme.It is recommended that you change Colab to the light theme.
The geemap package is pre-installed in Google Colab and is updated to the latest minor or major release every few weeks. Some optional dependencies of geemap being used by this notebook are not pre-installed in Colab. Uncomment the following line to install geemap and some optional dependencies.
# %pip install -U "geemap[workshop]"
Note that some geemap features do not work properly with Google Colab. If you are familiar with Anaconda or Miniconda, it is recommended to create a new conda environment to install geemap and its optional dependencies on your local computer.
conda create -n gee python=3.11
conda activate gee
conda install -c conda-forge mamba
mamba install -c conda-forge geemap pygis
Import the earthengine-api and geemap.
import ee
import geemap
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()
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
.
m = geemap.Map()
To display it in a Jupyter notebook, simply ask for the object representation:
m
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
.
Generate a map that focuses on the contiguous United States.
m = geemap.Map(center=[40, -100], zoom=4)
m
Generate a map that focuses on Alaksa.
m = geemap.Map(center=[64, -152], zoom=4)
m
Generate a map that focuses on University of Alaska Fairbanks.
m = geemap.Map(center=[64.864983, -147.840441], zoom=13)
m
To hide a control, set control_name
to False
, e.g., draw_ctrl=False
.
m = geemap.Map(data_ctrl=False, toolbar_ctrl=False, draw_ctrl=False)
m
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.
m = geemap.Map(basemap="Esri.WorldImagery")
m
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:
m.add_basemap("Esri.WorldTopoMap")
m.add_basemap("OpenTopoMap")
Print out the first 10 basemaps:
basemaps = list(geemap.basemaps.keys())
len(geemap.basemaps)
basemaps[:10]
You can also change basemaps interactively using the basemap GUI.
Google basemaps are not included in the geemap built-in basemaps due to license issues (source). Users can choose to add Google basemaps at their own risks as follows:
ROADMAP: https://mt1.google.com/vt/lyrs=m&x={x}&y={y}&z={z}
SATELLITE: https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}
TERRAIN: https://mt1.google.com/vt/lyrs=p&x={x}&y={y}&z={z}
HYBRID: https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}
m = geemap.Map()
url = "https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}"
m.add_tile_layer(url, name="Google Satellite", attribution="Google")
m
geemap.
in the code cell.
ee.
in the code cell.
m = geemap.Map()
m.add("basemap_selector")
m
Toggle the checkbox to show or hide the layer. Drag and move the slider to change the transparency level of the layer.
m = 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"],
}
m.add_layer(dem, vis_params, "SRTM DEM")
m.add_layer(states, {}, "US States")
m.add("layer_manager")
m
Click on the map to query Earth Engine data at a specific location.
m = geemap.Map(center=(40, -100), zoom=4)
dem = ee.Image("USGS/SRTMGL1_003")
landsat7 = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")
states = ee.FeatureCollection("TIGER/2018/States")
vis_params = {
"min": 0,
"max": 4000,
"palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}
m.add_layer(dem, vis_params, "SRTM DEM")
m.add_layer(
landsat7,
{"bands": ["B4", "B3", "B2"], "min": 20, "max": 200, "gamma": 2.0},
"Landsat 7",
)
m.add_layer(states, {}, "US States")
m.add("inspector")
m
m = geemap.Map(center=(40, -100), zoom=4)
dem = ee.Image("USGS/SRTMGL1_003")
vis_params = {
"min": 0,
"max": 4000,
"palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}
m.add_layer(dem, vis_params, "SRTM DEM")
m.add("layer_editor", layer_dict=m.ee_layers["SRTM DEM"])
m
m = geemap.Map(center=(40, -100), zoom=4)
landsat7 = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")
m.add_layer(
landsat7,
{"bands": ["B4", "B3", "B2"], "min": 20, "max": 200, "gamma": 2.0},
"Landsat 7",
)
m.add("layer_editor", layer_dict=m.ee_layers["Landsat 7"])
m
m = geemap.Map(center=(40, -100), zoom=4)
states = ee.FeatureCollection("TIGER/2018/States")
m.add_layer(states, {}, "US States")
m.add("layer_editor", layer_dict=m.ee_layers["US States"])
m
You can draw shapes on the map using the draw control. The drawn features will be automatically converted to Earth Engine objects, which can be accessed in one of the following ways:
ee.Geometry()
, use: m.user_roi
ee.FeatureCollection()
, use: m.user_rois
m = geemap.Map(center=(40, -100), zoom=4)
dem = ee.Image("USGS/SRTMGL1_003")
vis_params = {
"min": 0,
"max": 4000,
"palette": "terrain",
}
m.add_layer(dem, vis_params, "SRTM DEM")
m.add("layer_manager")
m
if m.user_roi is not None:
image = dem.clip(m.user_roi)
m.layers[1].visible = False
m.add_layer(image, vis_params, "Clipped DEM")
The Earth Engine Data Catalog hosts a variety of geospatial datasets. As of March 2024, the catalog contains over 1,100 datasets with a total size of over 100 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.
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. 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.
m = geemap.Map()
m
from geemap.datasets import DATA
m = geemap.Map()
dataset = ee.Image(DATA.USGS_GAP_AK_2001)
m.add_layer(dataset, {}, "GAP Alaska")
m.centerObject(dataset, zoom=4)
m
from geemap.datasets import get_metadata
get_metadata(DATA.USGS_GAP_AK_2001)
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.
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.
Images can be loaded by passing an Earth Engine asset ID into the ee.Image
constructor. You can find image IDs in the Earth Engine Data Catalog. For example, to load the NASA SRTM Digital Elevation you can use:
image = ee.Image("USGS/SRTMGL1_003")
image
m = 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"], # 'terrain'
}
m.add_layer(image, vis_params, "SRTM")
m
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.
For example, to load the image collection of the Sentinel-2 surface reflectance:
collection = ee.ImageCollection("COPERNICUS/S2_SR")
Let's find out how many Sentinel-2 images covering University of Alaska Fairbanks.
geometry = ee.Geometry.Point([-147.748432, 64.830691])
images = collection.filterBounds(geometry)
images.size()
images.first()
images = (
collection.filterBounds(geometry)
.filterDate("2023-07-01", "2023-09-01")
.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 5))
)
images.size()
m = geemap.Map()
image = images.first()
vis = {
"min": 0.0,
"max": 3000,
"bands": ["B4", "B3", "B2"],
}
m.add_layer(image, vis, "Sentinel-2")
m.centerObject(image, 8)
m
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:
m = geemap.Map()
image = images.median()
vis = {
"min": 0.0,
"max": 3000,
"bands": ["B8", "B4", "B3"],
}
m.add_layer(image, vis, "Sentinel-2")
m.centerObject(geometry, 8)
m
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.
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:
m = geemap.Map()
fc = ee.FeatureCollection("TIGER/2016/Roads")
m.set_center(-147.759590, 64.837553, 12)
m.add_layer(fc, {}, "Census roads")
m
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
feat = fc.first()
feat.toDictionary()
geemap.ee_to_df(fc)
m = geemap.Map()
states = ee.FeatureCollection("TIGER/2018/States")
fc = states.filter(ee.Filter.inList("NAME", ["California", "Oregon", "Washington"]))
m.add_layer(fc, {}, "West Coast")
m.center_object(fc, 5)
m
region = m.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)
m.add_layer(fc, {}, "Southeastern U.S.")
m.center_object(fc, 6)
m = geemap.Map(center=[40, -100], zoom=4)
states = ee.FeatureCollection("TIGER/2018/States")
m.add_layer(states, {}, "US States")
m
m = geemap.Map(center=[40, -100], zoom=4)
states = ee.FeatureCollection("TIGER/2018/States")
style = {"color": "0000ffff", "width": 2, "lineType": "solid", "fillColor": "FF000080"}
m.add_layer(states.style(**style), {}, "US States")
m
m = 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"]
m.add_styled_vector(
states, column="NAME", palette=palette, layer_name="Styled vector", **vis_params
)
m
Create a cloud-free imagery of Alaska. You can use Landsat 8/9 or Sentinel-2 imagery. Relevant Earth Engine assets:
m = 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}
m.add_layer(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,
}
m.add_layer(hyperion, hyperion_vis, "Hyperion")
m.add_plot_gui()
m
Set plotting options for Landsat.
m.set_plot_options(add_marker_cluster=True, overlay=True)
Set plotting options for Hyperion.
m.set_plot_options(add_marker_cluster=True, plot_type="bar")
from geemap.legends import builtin_legends
for legend in builtin_legends:
print(legend)
Add NLCD WMS layer and legend to the map.
m = geemap.Map(center=[40, -100], zoom=4)
m.add_basemap("Esri.WorldImagery")
m.add_basemap("NLCD 2021 CONUS Land Cover")
m.add_legend(builtin_legend="NLCD", max_width="100px", height="455px")
m
Add NLCD Earth Engine layer and legend to the map.
m = geemap.Map(center=[40, -100], zoom=4)
m.add_basemap("Esri.WorldImagery")
nlcd = ee.Image("USGS/NLCD_RELEASES/2021_REL/NLCD/2021")
landcover = nlcd.select("landcover")
m.add_layer(landcover, {}, "NLCD Land Cover 2021")
m.add_legend(
title="NLCD Land Cover Classification", builtin_legend="NLCD", height="455px"
)
m
Add a custom legend by specifying the colors and labels.
m = geemap.Map(add_google_map=False)
keys = ["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)]
m.add_legend(keys=keys, colors=colors, position="bottomright")
m
Add a custom legend by specifying a dictionary of colors and labels.
m = geemap.Map(center=[40, -100], zoom=4)
m.add_basemap("Esri.WorldImagery")
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/2021_REL/NLCD/2021")
landcover = nlcd.select("landcover")
m.add_layer(landcover, {}, "NLCD Land Cover 2021")
m.add_legend(title="NLCD Land Cover Classification", legend_dict=legend_dict)
m
Add a horizontal color bar.
m = geemap.Map()
dem = ee.Image("USGS/SRTMGL1_003")
vis_params = {
"min": 0,
"max": 4000,
"palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}
m.add_layer(dem, vis_params, "SRTM DEM")
m.add_colorbar(vis_params, label="Elevation (m)", layer_name="SRTM DEM")
m
Add a vertical color bar.
m.add_colorbar(
vis_params,
label="Elevation (m)",
layer_name="SRTM DEM",
orientation="vertical",
max_width="100px",
)
Make the color bar background transparent.
m.add_colorbar(
vis_params,
label="Elevation (m)",
layer_name="SRTM DEM",
orientation="vertical",
max_width="100px",
transparent_bg=True,
)
Create a split map with basemaps. Note that ipyleaflet has a bug with the SplitControl. You can't pan the map, which should be resolved in the next ipyleaflet release.
m = geemap.Map()
m.split_map(left_layer="Esri.WorldTopoMap", right_layer="OpenTopoMap")
m
Create a split map with Earth Engine layers.
m = geemap.Map(center=(40, -100), zoom=4, height=600)
nlcd_2001 = ee.Image("USGS/NLCD_RELEASES/2019_REL/NLCD/2001").select("landcover")
nlcd_2021 = ee.Image("USGS/NLCD_RELEASES/2021_REL/NLCD/2021").select("landcover")
left_layer = geemap.ee_tile_layer(nlcd_2001, {}, "NLCD 2001")
right_layer = geemap.ee_tile_layer(nlcd_2021, {}, "NLCD 2021")
m.split_map(left_layer, right_layer)
m
Create a 2x2 linked map for visualizing Sentinel-2 imagery with different band combinations. Note that this feature does not work properly with Colab. Panning one map would not pan other maps.
image = (
ee.ImageCollection("COPERNICUS/S2")
.filterDate("2023-07-01", "2023-09-01")
.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 5))
.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=[64.864983, -147.840441],
zoom=12,
ee_objects=[image],
vis_params=vis_params,
labels=labels,
label_position="topright",
)
Check the available years of NLCD.
m = 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. Note that ipyleaflet has a bug with the SplitControl. You can't pan the map, which should be resolved in the next ipyleaflet release.
m.ts_inspector(
left_ts=collection,
right_ts=collection,
left_names=years,
right_names=years,
left_vis=vis_params,
right_vis=vis_params,
width="80px",
)
m
Note that this feature may not work properly with Colab. Restart Colab runtime if the time slider does not work.
Create a map for visualizing MODIS vegetation data.
m = 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",
}
m.add_time_slider(collection, vis_params, time_interval=2)
m
Create a map for visualizing weather data.
m = 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)]
m.add_time_slider(collection, vis_params, labels=labels, time_interval=1, opacity=0.8)
m
Visualizing Sentinel-2 imagery
m = geemap.Map()
collection = (
ee.ImageCollection("COPERNICUS/S2_SR")
.filterBounds(ee.Geometry.Point([-147.840441, 64.864983]))
.filterMetadata("CLOUDY_PIXEL_PERCENTAGE", "less_than", 10)
.filter(ee.Filter.calendarRange(6, 8, "month"))
)
vis_params = {"min": 0, "max": 4000, "bands": ["B8", "B4", "B3"]}
m.add_time_slider(collection, vis_params)
m.set_center(-147.840441, 64.864983, 12)
m
Create a split map for visualizing the ESA WorldCover for the state of Alaska, with the Esri.WorldImagery
as the left layer, and the ESA WorldCover as the right layer. Add the ESA land cover legend to the map. Relevant Earth Engine assets: