Extracting pixel time series from Google Earth Engine

Question

I have written a short program in Python to extract a time series for any given pixel for MODIS data stored in the Google Earth Engine. The code is working fine and returns a data frame containing the relevant band value and date.

import pandas as pd
import numpy as np
from datetime import datetime as dt
import ee


def extract_time_series(lat, lon, start, end, product_name, band_name, sf):

    # Set up point geometry
    point = ee.Geometry.Point(lon, lat)

    # Obtain image collection for all images within query dates
    coll = ee.ImageCollection(product_name)\
        .filterDate(start, end)

    # Get list of images which correspond with the above
    images = [item.get('id') for item in coll.getInfo().get('features')]

    store = []
    date_store = []

    # Loop over all images and extract pixel value
    for image in images:

        im = ee.Image(image)

        # Obtain date from timestamp in metadata
        date = dt.fromtimestamp(im.get("system:time_start").getInfo() / 1000.)
        date_store.append(np.datetime64(date))

        # Extract pixel value
        data = im.select(band_name)\
            .reduceRegion(ee.Reducer.first(), point, 1)\
            .get(band_name)

        store.append(data.getInfo())

    # Scale the returned data based on scale factor
    store = [x * sf if isinstance(x, int) else np.nan for x in store]

    # Convert output into pandas data frame
    df = pd.DataFrame(index=date_store, data=store, columns=[band_name])

    return df


if __name__ == "__main__":

    ee.Initialize()

    latitude = 9.95
    longitude = -10.09 
    start_date = '2018-01-01'
    end_date = '2018-01-31'
    product = 'MODIS/006/MOD11A1' 
    band = 'LST_Day_1km'
    scale_factor = 0.02

    # Extract data and obtain pd.DataFrame
    output = extract_time_series(latitude,
                                 longitude,
                                 start_date,
                                 end_date,
                                 product,
                                 band,
                                 scale_factor)

However, I am aware of the way GEE handles reducing regions and how it can sometimes be difficult to know precisely what is happening. I decided to compare my extract MODIS values with the same time series requested from AppEEARS.

I noticed that the values are different and thus my code must be performing some spatial average, probably during the:

data = im.select(band_name)\
        .reduceRegion(ee.Reducer.first(), point, 1)\
        .get(band_name)

line. I set scale == 1 here so that the returned value I receive should be equal to the value at the native resolution (info on this).

How can I extract the MODIS band value for ONLY my point of interest?

Answer 1

Hope you find useful this tutorial:

http://www.loicdutrieux.net/landsat-extract-gee/examples.html

from geextract import ts_extract, get_date
from datetime import datetime
import numpy as np
import matplotlib.pyplot as plt
plt.figure(figsize=(10,5))

# Extract a Landsat 7 time-series for a 500m radius circular buffer around
# a location in Yucatan
lon = -89.8107197
lat = 20.4159611
raw_dict = ts_extract(lon=lon, lat=lat, sensor='LE7',
                      start=datetime(1999, 1, 1), radius=500)

# Function to compute ndvi from a dictionary of the list of dictionaries returned
# by ts_extract
def ndvi(x):
    try:
        return (x['B4'] - x['B3']) / (x['B4'] + x['B3'])
    except:
        pass

# Build x and y arrays and remove missing values
x = np.array([get_date(d['id']) for d in raw_dict])
y = np.array([ndvi(d) for d in raw_dict], dtype=np.float)
x = x[~np.isnan(y)]
y = y[~np.isnan(y)]

# Make plot
plt.plot_date(x, y, "--")
plt.plot_date(x, y)
plt.title("Landsat 7 NDVI time-series Uxmal")
plt.ylabel("NDVI (-)")
plt.grid(True)
plt.show()

Answer 2

I found the issue. It wasn't the scale, since anything below the native resolution of the product returns the value at the native resolution.

The problem was actually a missing parameter in the following line:

data = im.select(band_name)\
        .reduceRegion(ee.Reducer.first(), point, 1)\
        .get(band_name)

I changed this to:

        data = im.select(band_name)\
        .reduceRegion(ee.Reducer.first(),
                      point,
                      1,
                      crs=projection)\
        .get(band_name)

where projection is equal to the image projection (projection = im.projection().getInfo()['crs']**im.projection().getInfo()['crs']). By specifying the native resolution of the data product, the reduceRegion() method then returns the native value for my pixel location.

Answer 3

I think your problem problem lies in your scale, as you specifically say it should be 1M in your reduceRegion(). The nominal Scale for the MODIS/006/MOD11A1 is 1000M. Set scale to 1000 to see if it works. I can't test it out for you because I don't have the Python module for ee installed properly.