Retrieving and reprojecting MODIS satellite products

Data from the MODIS satellites are a great source for comparison to vegetation models and also for teaching. However, they come in a format I am not really used to and in a strange projection. The format is HDF, of which even several versions exists. But there are also several tools to convert it to a “better” file format for me. HDF needs much less space compared to other formats. Therefore it makes sense for satellite data to be stored as HDF.

First I need to know which data I want and in which region. Just recently a colleague asked me to get the fire product (MCD45A1) for Vietnam for him. Therefore I show it as example here.

On the map provided on the MODIS overview webpage one can try to select the correct row (v) and column (h) values. However, I already said it is a strange projection. I reprojected the grid to a standard equirectangular projection:

modis sinusoidal

That’s also not easy to select the right h/v values. However with both maps it worked. I set up a perl script using my module WebService::MODIS, which is also available on CPAN. Here is the example script, how to download the desired satellite data:

#!/usr/bin/env perl
use warnings;
use strict;
use WebService::MODIS;
readCache; # or "initCache;" if not done yet
my $fire = WebService::MODIS->new();
$fire->product("MCD45A1");
$fire->version('051');
$fire->h([27,28]);
$fire->v([6,7,8]);
$fire->dates(["2007.01.01", "2014.12.31"]);
$fire->createUrl;
#print "$_\n" foreach $fire->url;
$fire->download();

Now the files are in the current working directory. I do the mosaicing and reprojecting with ModisDownload.R. You may now ask: “Why not also downloading the data with the R Skript?” I did it in the beginning and it was very frustrating, since it crashed or did not start downloading at all and continued download was not possible. Therefore the perl module. To make use of different directories, where the downloaded files were and the reprojected files should be saved, I modified the original ModisDownload.R. Here is the diff of line 184 and 185 (Version 3.2, 13th Jan. 2014):

- write(paste(getwd(),"/",hdfNames[1], sep=""), mosaicname)
- for (j in 2:length(hdfNames)) write(paste(getwd(),"/",hdfNames[j], sep=""),mosaicname,append=T)
+ write(hdfNames[1], mosaicname)
+ for (j in 2:length(hdfNames)) write(hdfNames[j],mosaicname,append=T)

Now the filenames given to the mosaicHDF function must contain the full path to the files if called from another directory. ModisDownload.R uses the Modis Reprojection Tool (MRT), which is written in Java. Last year we used it for teaching. However, the installation under Windows it is a pain in the … Therefore, we preprocessed the data this year and supplied the students with GeoTiffs. Here is now the R script, which processes the downloaded hdf files and creates one GeoTiff per timeslice in UTM projection in the current working directory:

MRT.PATH <- "/path/to/MRT/bin"
MODISDOWNLOAD.DIR <- "/path/to/ModisDownload.R"
source(paste(MODISDOWNLOAD.DIR, "ModisDownload.R", sep="/"))
MODIS <- list(product="MCD45A1", bands="1 0 0 0 0 0 0 0", pixel.size=500)
DOWNLOAD.DIR <- "/path/to/downloaded/MODIS/files"
OUTPUT.DIR <- "/path/to/large/disk/for/GeoTIFFs"
setwd(OUTPUT.DIR)
downloaded <- list.files(DOWNLOAD.DIR, pattern="hdf$", full.names=TRUE)
if (length(downloaded) != 0) {
  yyyyddd <- sub("^.*\\.A", "", downloaded)
  yyyyddd <- unique(sub("\\..*$", "", yyyyddd))
  year <- trunc(as.integer(yyyyddd) / 1000)
  day.of.year <- as.integer(yyyyddd) - year * 1000
  dates <- as.Date(paste(year,"-01-01", sep=""), format="%Y-%m-%d") + day.of.year - 1
  for (j in 1:length(dates)) {
    in.hdf <- downloaded[grepl(paste("A", yyyyddd[j], sep=""), downloaded)]
    out.hdf <- paste(MODIS[["product"]], "_", dates[j], ".hdf", sep="")
    out.tif <- paste(MODIS[["product"]], "_", dates[j], ".tif", sep="")
    if (length(in.hdf)>1) {
      mosaicHDF(in.hdf, out.hdf, MRTpath=MRT.PATH)
    } else {
      file.symlink(in.hdf, out.hdf)
    }
    reprojectHDF(out.hdf,
      out.tif,
      MRTpath=MRT.PATH,
      proj_type="UTM",
      proj_params="0 0 0 0 0 0 0 0 0 0 0 0 0 0 0",
      utm_zone=0,
      datum="WGS84",
      bands_subset=MODIS[["bands"]],
      pixel_size=MODIS[["pixel.size"]])
    ## unlink(out.hdf)
  }
}

Now I have a timeseries of Geotiffs in UTM projection I can work with.

Advertisements

3 thoughts on “Retrieving and reprojecting MODIS satellite products

  1. Very nice scripts, thanks!
    Although I might ask why you don’t just use the hdf as provided? HDF5 (and also 4 to a lesser extent) are great formats which have excellent data compression, I/O speed and grouping capabilities.
    geoTiffs on the other hand…

    Accessing them in R is relatively straight forward using the rhdf5 package:
    https://www.bioconductor.org/packages/release/bioc/html/rhdf5.html
    Example:
    “`
    library(“rhdf5”)
    fname <- "yourdata.hdf"
    #List the content of the HDF5 file.
    tmp <- h5ls(fname)
    lon <- h5read(fname, "lon") # to get al the lon coordinates
    sst <- h5read(fname, "burndate")
    “`

    The MODIS package in R also has some functions to directly address hdf4 files
    https://conservationecology.wordpress.com/2014/08/11/bulk-downloading-and-analysing-modis-data-in-r/

    1. Hi Martin,

      I finally managed to import the MODIS hdf files directly into R. I didn’t want to install bioconductor and the hdf5 library doesn’t work, since the MODIS data is in HDF4. But in the end it was even easier than I thought: the raster package can import HDF4 directly. Here is an example:

      library(raster)
      data.dir <- /path/to/where/the/files/are/"
      fin <- c("MOD17A3.A2010001.h11v04.055.2011276120655.hdf",
      "MOD17A3.A2010001.h12v04.055.2011276121028.hdf")

      x <- raster(paste('HDF4_EOS:EOS_GRID:"', data.dir, fin[1], '":MOD_Grid_MOD17A3:Gpp_1km', sep=""))
      y <- raster(paste('HDF4_EOS:EOS_GRID:"', data.dir, fin[2], '":MOD_Grid_MOD17A3:Gpp_1km', sep=""))

      xy <- merge(x=x, y=y)
      spplot(xy)

      ll <- CRS("+proj=longlat +ellps=WGS84")
      z <- projectRaster(xy, res=30/3600, crs=ll, method="ngb")
      spplot(z)

      I only tried it on my Linux box and don't know if it will also be working in class where we are sadly working on Windows machines.

      Greetings
      Jörg

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s