2

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

Can anyone shed any light on how I can extract the MODIS band value for ONLY my point of interest? Any amendments to the code is also very welcome!

2

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

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.

  • Specifying any resolution lower than the native (i.e <1000m) doesn't make any difference. See here: developers.google.com/earth-engine/scale. The issue was that I needed to specify the crs for the reduceRegion() method. – tda Feb 22 at 14:29
1

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.

  • Strange that you need to specify the projection, as the docs specifically say crs (Projection, default: null): The projection to work in. If unspecified, the projection of the image's first band is used. If specified in addition to scale, rescaled to the specified scale. – Sean Roulet Feb 22 at 14:48
  • Yeah that's exactly what I thought... perhaps it's because the geometry is in lat/lon rather than the native projection. – tda Feb 23 at 10:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.