I have 11 files showing 8day daytime LST temps in .tif format from MODIS.

Is there any way to calculate the average 8day daytime temperature for march 2013 for example but using in the calculations only pixels which are not cloud-contaminated in any of the 11 files that participate in this process?


Is this the MOD11A2 product you're using? You have two issues to solve:

  1. How do you define "cloud-contamined" from the data?
  2. What toolset to process/average your scenes?

Let's say you have 11 files for the same grid tile, different 8-day periods, as GeoTIFF. For point 1. you could either use the Science Data Set layer "Clear_sky_days: the days in clear sky conditions and with valid LSTs" from your MOD11A2 product, or the MOD 35 cloud mask product. I'd use the first, as it is part of your dataset and produced for use with LSTs. Note that this layer is an 8-bit field with the ith bit marking the ith day as "clear" or "not clear". You only want days for which this value is 255 (11111111 in binary). Check out the LPDAAC website for further information.

Your basic process would be as follows: a) Load LST dataset and clear-sky days (looping through your files) b) construct a mask to mark all pixels that have any non-clear day, adding the pixels from your current dataset to this mask as you loop through them c) average your datasets only over the non-masked pixels

I would use Python with numpyand the GDAL bindings for this purpose. numpy has a masked array datastructure. You could progressively construct your masked array from the clear-sky days dataset:

from osgeo import gdal
current_clearsky = gdal.Open('current-clear-sky-tiff-file')
cclear_raster = current_clearsky.ReadAsArray()    
cclear_masked = ma.masked_less(cclear_raster, 255) # check actual values in clear-sky dataccset

You have various ways of how to assemble the data and the mask from here. See the documentation for numpy.ma and here. You could for example stack 11 masked arrays with the LST data in the data fields and the clear-sky masks as masks, and then average along the vertical dimension.


In addition to the smart and easy Python solution provided by @chryss, here is an identical approach using R along with the MODIS package. The package is currently being developed by Matteo Mattiuzzi and others and offers some very nice and useful utilities when working with MODIS data. See R-Forge for the latest package version and Matteo's ownCloud for the corresponding manual. There is also a tutorial by Steven Mosher published on R-bloggers that comes in handy when installing the package for the first time.

## Global settings

# Packages

# Working directory
path.wd <- "path/to/working/directory"

# MODIS paths
MODISoptions(localArcPath = paste(path.wd, "data/raw", sep = "/"), 
             outDirPath = paste(path.wd, "data/tif", sep = "/"))

I downloaded and processed MOD11A2 daytime LST for March 2013 and chose Egypt as geographic extent (for reasons of cloud cover minimization). However, you should bear in mind that 8-day composites, e.g. for March 2013, do hardly ever start or end exactly with the desired month, so a daily product (MOD11A1) might rather be what you are looking for.

## Data processing

# Download
time.range <- transDate(begin = "2013-02-26", end = "2013-03-31")
runGdal("MOD11A2", begin = time.range$beginDOY, end = time.range$endDOY, 
        extent = "Egypt", SDSstring = "100000000010", 
        job = "lst_egypt_mar13")

After downloading and extracting the data, I imported the single GeoTiffs using lapply along with raster. All daytime LST pixels that are not entirely cloud-free were then rejected based on the corresponding clear sky information as depicted by @chryss.

# Import
lst <- lapply(list.files("data/tif/lst_egypt_mar13", 
                         pattern = "LST_Day_1km", full.names = T), raster)
csd <- lapply(list.files("data/tif/lst_egypt_mar13", 
                         pattern = "Clear_sky_days", full.names = T), raster)

# Cloud-free pixels per 8-day interval
lst_cf <- foreach(i = lst, j = csd, .combine = "stack") %do% {
  i[j[] < 255] <- NA

Finally, I overlayed the quality-adjusted 8-day LST layers supplying mean as function argument, thus receiving (almost) monthly mean daytime LST. As you can see, only a small amount of the pixels covering the specified spatial extent are not contaminated by clouds in March 2013.

# Monthly average LST of cloud-free pixels
lst_cf_mth <- overlay(lst_cf, fun = mean)



  • Very nice, +1. Does the MODIS package also handle L1B and swath data? – chryss Sep 16 '13 at 18:50
  • Not yet, unfortunatelly. I recently contacted the development team because of that topic and they said it is planned, but not yet implemented. – fdetsch Sep 16 '13 at 20:13
  • Hm, that would be interesting. Thanks for the edit, btw. No idea why my fingers typed 288. – chryss Sep 16 '13 at 20:42

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.