Google Earth Engine use via Python, containers and other mythical beasts

The easiest way to run any Python stuff on a SE VM (actually an OSgeolive 15 image with some local stuff) is running Docker here. Which is also a modern way of running stuff in isolated environments aka containers.

Specifically, it allows to run a Jupytarlab instance, using an official image taken from Docker Hub.

All information are available here and the full documentation is https://jupyter-docker-stacks.readthedocs.io/en/latest/index.html.

Different flavors of Jupyter docker stacks can be installed and run under docker, they are illustrated here: https://jupyter-docker-stacks.readthedocs.io/en/latest/using/selecting.html

So one can add an official image to SE VM with:

docker pull jupyter/scipy-notebook

which is one of most complete Docker image and includes an Anacoda installation. After that one can use:

docker run -p 10000:8888 -v ${PWD}:/home/jovyan/work jupyter/scipy-notebook

note that the container run as non privileged user (ID=1001 in this case) all time, so you should take care of the ownership of the files.

First of all, note that the default jovyan user is not included in sudoers and has a disabled password. So, in order to run privileged commands one has to manage a few chages by using a root interactive session on the running container:

docker ps
docker exec -it -u 0 <container-id> /bin/bash
[ ... make changes ... ]
docker commit <container-id> <new-image-name>

Within this session it is possible to add a few useful stuff to the base image, e.g.

apt install gdal-bin gdal-dev
python3 -m pip install rasterio

Let’s try to use Google Earth Engine directly in the SE VM, it could conflict with other stuff, but we could eventually revert to use Docker in case of problems.

The basic use of Python API is very similar to the Javascript one.

#!python3 -m pip install earthengine-api

#!python3 -m pip install earthengine-api --upgrade

import ee
ee.Authenticate()
ee.Initialize()

To authorize access needed by Earth Engine, open the following URL in a web browser and follow the instructions:

https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=tYrbiziekwJR5KGf8BSkdFz7V1VPUdU_YQn23qKzNgY&tc=jtvL0uK4Ko5K50yD7mYdSYSjC-R7v76VnmCvlnBmhx0&cc=POqtZ-KRzgCTNgvrBUTxf-DVnUmy1nGid2GD16c2izg

The authorization workflow will generate a code, which you should paste in the box below.

Enter verification code: 4/1AbUR2VMAhSyJisx3XRgzXUIu5RHLkY1xnwsIWJOr-Agl2jmeK1J_neggo40

Successfully saved authorization token.

import folium

lat, lon =  45.77, 4.855

mymap = folium.Map(location=[lat, lon], zoom_start=10)
mymap

Make this Notebook Trusted to load map: File -> Trust Notebook

An easier way of integrating graphics to the regular ee class methods is using the Geemap third-party package. It can use both folium or leaflet Javascript clients. The choice of the client depends on the status of their support.

#!python3 -m pip install geemap

#!python3 -m pip install jupyter_contrib_nbextensions

import geemap

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

You need to autheticate (eventually with 2FA) every time you use GEE, even via Python. Note that a one time token is typically used, but it is required to authorize the Jupyter client for use of GEE.

If the ipyleaflet widget does not appear, typically you need to restart the IPython kernel.

Map.add_basemap('Esri.OceanBasemap')

Map.add_basemap('Esri.NatGeoWorldMap')

Map

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

Map2

This is the common way to use folium instead of leaflet in the notebook. Currently seems the most stable way of using Geemap.

import geemap.foliumap as geemap
Map = geemap.Map(center=[40,-100], zoom=4)

Map

Map2 = geemap.Map(center=[40,-100], zoom=5)
Map2.add_basemap('Esri.OceanBasemap')
Map2.add_basemap('Esri.NatGeoWorldMap')

Map2

Display an NDWI image

import folium

image = ee.Image('LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318')

# Create an NDWI image, define visualization parameters and display.
ndwi = image.normalizedDifference(['B3', 'B5'])
ndwi_viz = {'min': 0.5, 'max': 1, 'palette': ['00FFFF', '0000FF']}

# Define a map centered on San Francisco Bay.
map_ndwi = folium.Map(location=[37.5010, -122.1899], zoom_start=10)

def add_ee_layer(self, ee_image_object, vis_params, name):
    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)
    folium.raster_layers.TileLayer(
      tiles=map_id_dict['tile_fetcher'].url_format,
      attr='Map Data &copy; <a href="https://earthengine.google.co/">Google Earth Engine</a>',
      name=name,
      overlay=True,
      control=True
  ).add_to(self)

folium.Map.add_ee_layer = add_ee_layer


# Add the image layer to the map and display it.
map_ndwi.add_ee_layer(ndwi, ndwi_viz, 'NDWI')
display(map_ndwi)

Make this Notebook Trusted to load map: File -> Trust Notebook

Zonal statistics

import ee
import geemap
import os

Ma = geemap.Map()
Map

dem = ee.Image('USGS/SRTMGL1_003')

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

out_dir = os.path.join(os.path.expanduser('~'), 'Downloads')
out_dem_stats = os.path.join(out_dir, 'dem_stats.csv')

if not os.path.exists(out_dir):
    os.makedirs(out_dir)

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

out_landsat_stats = os.path.join(out_dir, 'landsat_stats.csv')
geemap.zonal_statistics(
    landsat, states, out_landsat_stats, statistics_type='SUM', scale=1000
)

geemap.create_download_link(out_dem_stats)

geemap.create_download_link(out_landsat_stats)

Convert JS to Python

import ee
import geemap

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

js_snippet = """

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

// Create an NDWI image, define visualization parameters and display.
var ndwi = image.normalizedDifference(['B3', 'B5']);
var ndwiViz = {min: 0.5, max: 1, palette: ['00FFFF', '0000FF']};
Map.addLayer(ndwi, ndwiViz, 'NDWI', false);

"""

geemap.js_snippet_to_py(js_snippet)

import ee

# Load an image.
image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318')

# Create an NDWI image, define visualization parameters and display.
ndwi = image.normalizedDifference(['B3', 'B5'])
ndwiViz = {'min': 0.5, 'max': 1, 'palette': ['00FFFF', '0000FF']}
Map.addLayer(ndwi, ndwiViz, 'NDWI', False)

Map

Map

Importing stuff

import ee
import geemap

geemap.ee_search()