11

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?

Improve this question
0

3 Answers 3

Reset to default
6

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()
Improve this answer
4

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.

Improve this answer
2
  • 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. Commented Feb 22, 2019 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. Commented Feb 23, 2019 at 10:14
2

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.

Improve this answer
1
  • 1
    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. Commented Feb 22, 2019 at 14:29

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.