Getting Started
Setting up accounts and testing some python.
This tutorial is all about getting everything setup so you can do the other tutorials on this site. In the process we will be reviewing fundamental python skills required. If at any point you're not sure what is happening please have a look at our background pages on Python, Geospatial Python, and Machine Learning.
Accounts
For the lessons on this site you will need a Google Account (this can be @gmail or an institution/custom google account). You will use this account to access Google Drive, Google Colab, and activate Google Earth Engine.
Google Earth Engine
You will also need an account for Google Earth Engine for some lessons. You can sign up here, using your Google Account from the previous step to link everything together nicely.
Zindi
For some lessons we use data from Zindi.africa a Machine Learning competition site. Sign up here.
from google.colab import drive
drive.mount('/content/drive/')
Explore your drive
Now that you've connected (mounted) google drive to your session, you can access it as a local disk.
Look at the panel on the left, for the folder icon (3rd from the top). You should see your google drive under drive where you should see My Drive and Shared Drives. At any point if you hover over a folder three vertical dots appear on the right, and if you click you will see the option to copy path.
Let's now explore using python to list files.
import os
# list the folders and files in your drive
mydrive = '/content/drive/My Drive'
colab = '/content/drive/My Drive/Colab Notebooks'
print(os.listdir(mydrive))
# list the current working directory
print(os.getcwd())
# make a data folder in your Colab Notebook directory if it doesn't exist
data_folder = os.path.join(colab,'data')
os.path.isdir(data_folder)
if (not os.path.isdir(data_folder)):
os.mkdir(data_folder)
# list the contents
print(os.listdir(colab))
Google Earth Engine
Initialize
Next you will need to install the Google Earth Engine python API. If running on Colab it's already installed. It will autodetect and skip if already installed.
Then you need to initialize the Earth Engine session. The output will show a link you have to open, copy the code from the page that loads and paste back into a cell provided.
# If not on Colab you'll need install the earth-engine Python API
#!pip install earthengine-api #earth-engine Python API
# Athenticate to your GEE account.
!earthengine authenticate
# Earth Engine Python API
import ee
ee.Initialize()
Define a Bounding Box
You can define a Rectangle, Point, or Polygon for a spatial search by typing numbers, but that's tedious. It's much easier to get it from an existing data set or by selecting on a map. Let's get a Bounding Box in geojson format from geojson.io
- Search for a place (eyeglass in upper right), e.g Kathmandu
- Copy and Paste the geojson from the right panel, or Save GeoJSON file and then upload it to your Google Drive.
Note: Due to some limitations we can’t use ipyleaflet in Colab so we needed to use an external tool.
# Make sure we have the libraries we need
%pip install folium
%pip install geopandas
# Let see it on a map
import folium
import geopandas as gpd
# Paste the Geojson from geojsonio
aoi_geojson = '''{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
85.26712417602539,
27.66558933380944
],
[
85.38213729858398,
27.66558933380944
],
[
85.38213729858398,
27.747202035778272
],
[
85.26712417602539,
27.747202035778272
],
[
85.26712417602539,
27.66558933380944
]
]
]
}
}
]
}'''
# Now let's read the geojson into a GeoPandasDataframe
aoi = gpd.read_file(aoi_geojson)
print(aoi) # so we can see what it looks like
#Get the bounding box
bbox = aoi.total_bounds
print(bbox) # see the coordinates
#Make it a GEE rectangle
gee_aoi = ee.Geometry.Rectangle(bbox.tolist())
center = aoi.centroid[0]
m = folium.Map(location=[center.y, center.x], tiles="OpenStreetMap", zoom_start=12)
folium.Marker(
location=[center.y, center.x],
popup='Center',
icon=folium.Icon(color='green', icon='ok-sign'),
).add_to(m)
folium.features.GeoJson(aoi_geojson,
style_function = lambda x: {'color':'blue', 'fillcolor':'transparent'}
).add_to(m)
m
Query Google Earth Engine
Now lets pick a collection of imagery to query for the region. We can choose from anything in the Google Earth Engine data catalog
Let's start with the Landsat 8 Surface Reflectance this year, using a cloud mask and selecting the bands to make an R,G,B mosiac
# To make a map we first need some helper functions
# Define the URL format used for Earth Engine generated map tiles.
EE_TILES = 'https://earthengine.googleapis.com/map/{mapid}/{{z}}/{{x}}/{{y}}?token={token}'
#@title Mapdisplay: Display GEE objects using folium.
def Mapdisplay(center, dicc, Tiles="OpensTreetMap",zoom_start=10):
'''
:param center: Center of the map (Latitude and Longitude).
:param dicc: Earth Engine Geometries or Tiles dictionary
:param Tiles: Mapbox Bright,Mapbox Control Room,Stamen Terrain,Stamen Toner,stamenwatercolor,cartodbpositron.
:zoom_start: Initial zoom level for the map.
:return: A folium.Map object.
'''
mapViz = folium.Map(location=center,tiles=Tiles, zoom_start=zoom_start)
for k,v in dicc.items():
if ee.image.Image in [type(x) for x in v.values()]:
folium.TileLayer(
tiles = v["tile_fetcher"].url_format,
attr = 'Google Earth Engine',
overlay =True,
name = k
).add_to(mapViz)
else:
folium.GeoJson(
data = v,
name = k
).add_to(mapViz)
mapViz.add_child(folium.LayerControl())
return mapViz
# Query GEE for Landsat
l8_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')\
.filterDate('2020-01-01', '2020-08-31')\
.median()
# Time to make a map
l8_vis_params = {
'bands': ['B4', 'B3', 'B2'],
'min': 0,
'max': 3000,
}
Mapdisplay(center=[center.y, center.x],
dicc={'L8':l8_image.getMapId(l8_vis_params)},
zoom_start=12)
Exporting
At some point we often want to export data from Earth Engine for other uses. One might ask, can't we just do the whole analysis in Earth Engine? In some cases yes that's possible, but in other cases you need access to algorithms, data, or map making tools that just aren't possible in Earth Engine.
# But exporting the whole world to work with would be a pain, lets subset to our aoi
band_sel = ('B2', 'B3', 'B4', 'B5')
l8_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')\
.filterDate('2020-01-01', '2020-08-31')\
.select(band_sel)\
.filterBounds(gee_aoi)\
.median()
l8_image.getInfo()
# Only run this when necessary, once exported you can use the file from google drive.
# Landsat 8 in this example takes about 10 minutes or less.
task.start()
%pip install rasterio
#%pip install geopandas # remember we already did this earlier
import rasterio as rio
import geopandas
%matplotlib inline
# Loading a raster and making a map
raster_path = os.path.join(output_mount_folder, output_drive_folder, raster_name)
l8_raster_path = ".".join([raster_path, "tif"])
# Check some attributes
with rio.open(l8_raster_path, 'r') as l8_raster:
print(l8_raster.bounds)
print(l8_raster.meta)
import matplotlib.pyplot as plt
from rasterio.plot import show
%matplotlib inline
with rio.open(l8_raster_path, 'r') as l8_raster:
fig, (ax1, ax2, ax3, ax4) = plt.subplots(ncols=4, nrows=1, figsize=(10, 4), sharey=True)
# Plot Red, Green and Blue (rgb)
show((l8_raster, 1), cmap='Blues', ax=ax1)
show((l8_raster, 3), cmap='Reds', ax=ax3)
show((l8_raster, 2), cmap='Greens', ax=ax2)
show((l8_raster, 4), cmap='magma', ax=ax4)
# Add titles
ax1.set_title("Blue")
ax2.set_title("Green")
ax3.set_title("Red")
ax4.set_title("NIR")
from rasterio.plot import reshape_as_image
import numpy
# Function to normalize the grid values
def normalize(array):
"""Normalizes numpy arrays into scale 0.0 - 1.0"""
array_min, array_max = array.min(), array.max()
return ((array - array_min)/(array_max - array_min))
# Open and read the raster data
with rio.open(l8_raster_path, 'r') as l8_raster:
l8_data = l8_raster.read()
# To plot as RGB we have to normalize the data
l8_image = numpy.empty(l8_data.shape, dtype=numpy.float)
for band in range(l8_data.shape[0]):
l8_image[band] = normalize(l8_data[band])
fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(10, 4), sharey=True)
show(l8_image[[2,1,0],:,:], transform=l8_raster.transform, adjust='linear', ax=ax1) # RGB
show(l8_image[[3,2,1],:,:], transform=l8_raster.transform, adjust='linear', ax=ax2) # False Color IR
