Tutorial 11: Python and Google Earth Engine

Some notes on documentation

In the documentation we will refer to today, you will notice that there are examples available for both JavaScript and Python (Google Colab). This is because Google is working on implementing all the JavaScript functionality in Python. Not all the documentation or examples are at this time available in Python. Luckily, although the languages are quite different, the JavaScript and Python implementation for this package are very similair. It is a good skill to be able to read some JavaScript. Additionally, we are not sure what documentation is available for python at the time of following this tutorial. Therefor we will refer to the JavaScript documentation, but we will highlight any important differences and provide guidance on using this code in Python.

Some common differences are the following. To initialize a variable, in Python you can simply type:

variablename = 'value'

In Javascript we can initialize a variable in various ways. The difference is where this variable in the code will be known:

var variablename = 'value'
let variablename = 'value'
const variablename = 'value'

Another often occuring difference in the documentation is how variables are passed to a function. In python we can do that as we are used to:

multiStats = img.reduceRegions(
  collection=regionCol,
  reducer=reducer,
  scale=30,
  crs='EPSG:3310',
)

In Javascript we have to pass an object, similar to a dictionary in python:

var multiStats = img.reduceRegions({
  collection: regionCol,
  reducer: reducer,
  scale: 30,
  crs: 'EPSG:3310',
});

You do not need to understand this code itself for now, just note the differences in syntax.

Initialize the map

By following these steps, we will create a map visualization of the selected image.

# set our initial zoom position for the netherlands
center_lon = 51.962589
center_lat = 5.669627
zoomlevel = 8

# The following line will initialize the map
Map = geemap.Map(center=[center_lon, center_lat], zoom=zoomlevel)

# The following line will actually call and show the map:
Map

The map will stay and remain in this window unless we instruct it otherwise. It is an interactive map with 1 standard basemap at the moment. To start over, you can call Map.clear(), which will remove all layers and other objects from the map, requiring you to initialize it again. Sometimes, when working with an interactive widget like this map, Jupyter might become unresponsive. If you no longer receive any response (you can test this by running 1 + 1 in a new code block), you might need to restart the kernel and begin again.

Select 1 image and add it to the map

Step 1 is complete! The next step is to select one image from the image collection and add it to the map.

# Initial date of interest (inclusive).
from_date = '2017-01-01'
# Final date of interest (exclusive).
to_date = '2023-01-01'

image = viirs_image_collection\
            .filterDate(from_date, to_date)\
            .sort("system:time_start", True)\
            .first()

visualization_params = {
    'min': 0,
    'max': 10,
    'opacity': 0.5,
    'bands': ['avg_rad']
}

Map.addLayer(ee_object = image, 
             vis_params = visualization_params,
             name = 'My First Image')

In the code snippet provided, several steps are performed to add an image to the map. In the first line the following happened:

1. Filtering Images: The viirs_image_collection is filtered to select only the images between two specified dates using the filterDate(from_date, to_date) method.

2. Sorting the Collection: The filtered image collection is sorted based on the "system:time_start" property in ascending order using the sort("system:time_start", ascending) method. The resulting collection is then assigned to the variable image.

3. Selecting an image: We selected the first image out of the collection using the first() method.

We then define visualization parameters and add the image to the map:

1. Visualization Parameters: Visualization parameters are defined to specify how the image should be displayed on the map. This includes the desired range of values, opacity, and the specific band to visualize (avg_rad in this case).

2. Adding the Image to the Map: The image is added to the map using the AddLayer method. This method requires the Earth Engine object (image), the visualization parameters, and an optional name for the layer.

By following these steps, the code successfully adds the selected image to the map, great!

Add the boundaries of the Netherlands to the map

The next code will add the boundaries of the Netherlands to the map. To get the boundary data we use the pygadm package. pygadm is a dataset project from the University of Berkeley that provides administrative boundary data on multiple scales for various regions worldwide. Be careful! The returned data is not always completely correct, but it suits our purposes for now.

# Get data
netherlands = pygadm.get_items(admin='NLD', content_level=1)

# Define the projection
netherlands.crs = "EPSG:4326"

# Import gpd dataframe to earth engine
ee_netherlands = geemap.geopandas_to_ee(netherlands)

# Add data to map
Map.addLayer(ee_object = ee_netherlands, 
             name = 'Netherlands Boundaries')
Map.addLayerControl()

Using the get_items function, we can obtain the boundaries at content level 1, which represent the provinces in the Netherlands. This data is retrieved in the GeoPandas GeoDataFrame format, which cannot be added to the map directly. However, we can utilize a convenient conversion tool provided by the geemap package to convert it into an Earth Engine feature. The resulting object can then be easily placed on the map.

For additional examples and tools, we recommend checking out the examples provided by geemap themselves! It’s worth noting that you can also export a map to HTML, which might be useful for your project.

Take a look at these examples to further enhance your project and explore the capabilities offered by geemap.

Great, we have now made a simple map in google earth engine. We added an image for the whole world and we imported external data and also visualized that. And it works very quick!

Selecting the data and collapsing it

The first step can be done in one line. Several times we call a method on an image collection, which will result in an image collection. On this collection we will call another method and so fort. Look at the following line of code and try to understand what is happening before you read the explanation. Use the documentation!

viirs_image = viirs_image_collection.select('avg_rad') \
                .filterDate(from_date, to_date) \
                .toBands() 

4 methods are called in one line. In this code we go from an image collection over all available time steps to a single image of the time we are interested in, with one band representing an image per month. Why we collapse it into this single image will become clear in the next step. The separate methods are described below.

1. First, we select the band we are interested in, the avg_rad band, containing the average nighttime radiation per month. For this we use the select method.
2. We then select the images from the dates we are interested in using the filterDate method.
3. Finally we call the toBands method. Note that this method results in an Image instead of an ImageCollection.

When you run the code above note that the code is run instantly, no matter the processing power of your laptop or virtual machine, while actually it is quite a large computation. How can this be? Actually the code and the calculation is not done on your computer, in fact there is no computing done at all yet for this code. This code creates an API Call, which will be sent to the Earth Engine cloud platform once it needs to be called. For example when you now call the viirs_image variable, it takes longer. The API call is sent and the results are returned to your environment.

viirs_image

Clipping the data

Next, we will clip the image to the desired area, the boundary of the Netherlands. This can be done by using the clip method. We could have also added this to the already long line of code above. Functionally, this makes no difference. The decision on how to do this is personal; every programmer might do it differently. The most important aspect of this choice is to make it clear for yourself and others what each line is doing and to define clear variable names.

# Function to clip an image to the specified geometry
viirs_image_clipped = viirs_image.clip(ee_netherlands)

viirs_image_clipped now contains an image with 72 bands. Please verify this in your own environment.

Calculating statistics

Now, we want to calculate the average radiance per province. We can achieve this using a reducer, which we apply over an area. A reducer is an algorithm that is applied over a set of values. In our case, the set of values will be the pixel values within a province. We will use a method called reduceRegions to gather this set of values. According to the documentation, this method requires a collection and a reducer as inputs. We will provide the polygons of the provinces in the Netherlands, referred to as ee_netherlands, as the input collection. As for the reducer, there are various options available. To explore the available options, you can search for reducer using the filter in the documentation.

Since we want to take the average over the collection, we will utilize the ee.Reducer.mean() function. It is important to note that this function itself is not a reducer, but rather it returns a reducer. Therefore, we need to actually call the function as part of the input. Additionally, we provide the optional parameter scale with a value of 463.83 meters.

province_statistics = viirs_image_clipped.reduceRegions(
    collection = ee_netherlands,
    reducer = ee.Reducer.mean(),
    scale = 463.83
)

Exporting

Geemap provides useful tools for exporting data. The most convenient format for us is a GeoPandas object, but there are many other options available as well. Please refer to the documentation for additional options.

df = geemap.ee_to_geopandas(province_statistics)

It is, obviously, also possible to export other earth engine objects like images or image collections, as is shown in the documentation. At the end of this tutorial we will show how to export the data used in this tutorial to a geotiff file.

Calculating differences

From here calculating the difference is straight forward, since we know (geo)pandas well. We now also have all the functionality that comes with GeoPandas, such as very easily writing the results to a geojson!

df['diff'] = df['20170101_avg_rad'] - df['20221201_avg_rad']
df.to_file('test.geojson', drive='GeoJSON')

Classifying and visualizing trends

Now to see how the radiation behaved over time we have to do a little more work, since the dataframe in this format is not very useful for this. Read and understand the code below.

import pandas as pd
import numpy as np 
import matplotlib.pyplot as plt
import matplotlib

# Clean the dataframe a little
clean_df = df.drop(['GID_0', 'GID_1', 'ISO_1', 'HASC_1', 'NAME_0', 'NL_NAME_1', 'VARNAME_1', 'TYPE_1','CC_1', 'ENGTYPE_1'], axis=1)

# Reformat the dataframe so the rows depict measurements per timestamp per province, instead of an attribute per timestamp
melted_df = pd.melt(clean_df, id_vars=['NAME_1', 'geometry', 'diff'] , var_name='date', value_name='value')

# Convert the time string to a datetime format
melted_df['short_date'] = melted_df['date'].str[:8].str.strip()
melted_df['datetime'] = pd.to_datetime(melted_df['short_date'], format="%Y%m%d")


# Initialize the plot
fig, ax = plt.subplots(figsize=(9, 6))

# Create an empty list to save all lines in for the legend 
lines = []

# For each province, add a line with a different color 
for prov in melted_df['NAME_1'].unique():    
    
    # Create a dataframe for the province of interest
    d = melted_df[melted_df['NAME_1'] == prov]
    
    # Generate a random color
    random_rgb_values = list(np.random.random(3))
    random_color = matplotlib.colors.rgb2hex(random_rgb_values, keep_alpha=False)
    
    # Add a line for each province
    lines += ax.plot( d["datetime"].values, d["value"].values, color= random_color, lw=1, zorder=10, label=prov)
    
    
# Add the legend    
ax.legend(handles = lines)

We can see that there are some differences per province, but that there does not seem to a nationwide trend. We can also see that in the summer there is no data for the Netherlands, or in fact for entire northern europe. This has to do with the nights being so long here in the summer. Per province there might be a larger difference. To find out, we can add a trendline to the graph, but this is out of scope for this tutorial.

Classify differences

For the visualization we now want to categorize the differences in 3 classes, a little difference, medium difference and a large difference. To do this we will use a numpy function called linspace and a pandas function called cut. This results in an array with datatype CategorizalClasses, which earth engine does not understand so we convert it to a string.

# Create 3 classes
breaks = np.linspace(df['diff'].min() - 0.0001, df['diff'].max(), 4)

# Assign values to classes
classes = ['Low', 'Mid', 'High']
df['class'] = pd.cut(df['diff'], bins=breaks, labels=classes, right=True)

# Convert the datatype to strings so earth engine understands
df['class'] = df['class'].astype('str')

Show results on a map

So the next step is to show the differences on an earth engine map. To do this, we have to convert the geopandas dataframe back to a earth engine object. Before we do this we have to tell geopandas what coordinate refrence system the geometries are in. Strangely the coordinate system is not exported from earth engine, but it does need it to be assigned to be imported…

df = df.set_crs(4326)
new_ee_statistics = geemap.geopandas_to_ee(df)

# If you run the entire tutorial in one go, the map will show at the same location as it was first initialized.
# To show it hear, clear the map and create a new one. 
# If the map was not initizalized before, this line will throw an error. 
Map.clear()

Style the map

Now we have the polygons back in the right format with the classes. To visualize we define a palette, which tells earth engine which color we will give to each class. For the legend we can use this dictionary directly. Of course we also want to add the latest night time imagery behind it and we want to add a layer control.

Map = geemap.Map(zoomlevel=8)

palette = {
    'High': 'B81D13', # red
    'Mid': 'EFB700',  # yellow
    'Low': '008450'   # green
}


visualization_params = {
    'min': 0,
    'max': 10,
    'opacity': 0.7,
    'bands': ['20221201_avg_rad']
}

Map.addLayer(ee_object = viirs_image_clipped, 
             vis_params = visualization_params,
             name = 'Night Time')

Map.add_styled_vector(new_ee_statistics, column="class", palette=palette, layer_name="Styled vector")
Map.center_object(new_ee_statistics)

Map.add_legend(legend_title="Legend", legend_dict=palette, position = 'topright')

Map.addLayerControl()

Map

To file

As promised earlier in this tutorial we will show how to export raster data to a local file for visualization or analysis outside of GEE. We will use the tools as described in the geemap tutorials. Before we export an image we have to get some data from GEE. We will use the same ata as used in the rest of the tutorial.

# Initial date of interest (inclusive).
from_date = '2017-01-01'
# Final date of interest (exclusive).
to_date = '2023-01-01'

image = viirs_image_collection\
            .filterDate(from_date, to_date)\
            .sort("system:time_start", True)\
            .first() 

# Get geometry data and convert to earth engine object
netherlands = pygadm.get_items(admin='NLD', content_level=1)
netherlands.crs = "EPSG:4326"
ee_netherlands = geemap.geopandas_to_ee(netherlands)

After collecting the data we can use a simple command for downloading the data. In this case the image we are working with consists of two bands. Setting the file_per_band parameter to True will result in a file for each of the bands.

# Export to geotiff
geemap.ee_export_image(image, 
                       filename='clipped_image.tif', 
                       region=ee_netherlands.geometry(), 
                       scale=285,
                       file_per_band=True
                       )

To python

For further analysis within python it might be easier to convert a raster to a numpy.array. This stackexchange question and answer by Justin Braaten shows nicely how this can be done. Notice that for further spatial analysis this might be inconvenient, since the resulting numpy array contains only the pixel values, without any spatial information. When importing to, for example, rasterio we need additional sptial information such as the origin (north west corner) and the pixel size, however. For more information see also this stackexchange thread