Download and resampling of MODIS images

From spatial-analyst.net
Jump to: navigation, search
T

his article explains how to automate download, mosaicking, resampling and import of MODIS product to a GIS. We focus on the one of the most known MODIS products for terrestrial environmental applications: the Enhanced Vegetation Index (EVI), which is the improved NDVI (Huete et al., 2002; see the complete list of MODIS products). EVI corrects distortions in the reflected light caused by the particles in the air as well as the ground cover below the vegetation. The EVI also does not become saturated as easily as the NDVI when viewing rainforests and other areas with large amounts of chlorophyll. EVI can be directly related to the photosynthetic production of plants, and indirectly to the green biomass (Huete et al., 2002). By observing dynamics of EVI for an area of Earth's surface, we can infer about the vegetation dynamics within season, but also detect long-term trends and sudden changes in the biomass (e.g. due to forest fires, deforestation, urban growth and similar).

NASA's MODIS Earth observation system is today considered to be one of the richest sources of remote-sensing data for monitoring of environmental dynamics (Neteler, 2005; Lunetta et al. 2006; Ozdogana and Gutman 2008). This is due to the following reasons:

  1. it is global;
  2. it has a relatively high temporal resolution/coverage (1-2 days);
  3. it is open-access (see the MODIS licence specification);
  4. a significant work is done to filter the the original raw images for clouds and artifacts; a variety of complete MODIS products such as composite 15-day and monthly images is available at three resolutions: 250 m, 500 m and 1 km;
  5. efficient tools to obtain various MODIS products and import them to various GIS exist;

The importance of MODIS EVI index for environmental management and nature conservation is enormous (Lunetta et al. 2006). MODIS images can be obtained by anybody at any time and can be used as an independent source of information to prove degradation of natural systems that are possibly not visible from local sources of information. MODIS EVI images can be used to prove changes in land cover, deforestation, damages caused by global warming or even fine long-term succession processes. Because MODIS EVI is available from year 2000, anybody can compare the current situation with the situation from a decade ago.

MODIS download and use via R.
  • Icon R.png MODIS_download.R : R script to automatically download, mosaick and resample MODIS blocks to some local projection system.
  • Ge icon.jpg MODIS_tiles.kmz : MODIS tiling system world map (based on the Gpolygon field in the XML files attached to the original tiles).
  • Icon zip.png MRT_download_Win.zip : MODIS resampling tool.


By adjusting this R script, you should be able to automate download and preparation of any type of MODIS products for your own mapping application. If you experience problems, please report them using the R-sig-geo mailing list, or start a discussion via this wiki.


Contents

Installation of software

Before you can start downloading MODIS EVI images from R, you need to obtain and install some necessary packages and R libraries:

  1. RCurl --- allows you to list directories on an FTP server (read more);
  2. wget --- this will automate download of images from R; this exe you only need to put in your Windows directory folder (Note: make sure you disable your antivirus tools such as Norton or McAfee otherwise it might block wget from running);
  3. MRT --- MODIS reproject tool can be used to mosaic MODIS images, resample them to other coordinate systems, and export images to more common GIS formats;

After you have finished installing these two packages, you also need to specify the location of MRT and directory where you want to output all EVI maps you will produce:

 
 > library(RCurl) 
 > library(rgdal)
 > MOD13A3 <- "ftp://e4ftl01.cr.usgs.gov/MOLT/MOD13A3.005/"
 > MOD13A3a <- "ftp://anonymous:test@e4ftl01.cr.usgs.gov/MOLT/MOD13A3.005/"
 > MRT <- 'E:\\MODIS\\MRT\\bin\\'  # MODIS Reprojection Tool exe files
 > workd <- 'E:\\PineBeetleBC\\MODIS\\'   # working directory
 

Download of MODIS HDF tiles

MODIS images are typically distributed as HDF (Hierarchical Data Format) 10 by 10 arcdegree-tiles, projected in the Sinusoidal projection (Sinusoidal projection has been promoted by geographers as the most suited projection for global image databases; see e.g. Chang Seong et al., 2002), both HDF format and Sinusoidal projection are not supported in many GIS. Therefore, before you can use MODIS images, you will need to run some pre-processing to glue the tiles and coerce the data in a more usable format. We are interested to obtain monthly 1 km resolution MODIS EVI images (this dataset is known by the name "MOD13A3"), resample them to our local projection system (EPSG:3005), and then export them to a more common GIS format, namely GeoTiff.

The MOD13A3 tiles can be browsed/obtained directly via NASA's FTP server. We start by fetching the list of sub-directories (the dates) of interest:

 
# get the list of directories (thanks to Barry Rowlingson):
 > items <- strsplit(getURL(MOD13A3), "\n")[[1]]
 > items[2]
 [1] "drwxr-xr-x  2 90 118784 Jan  5  2009 2000.02.01\r"
# you get the folders (and files) but the folder names are in the form of a unix directory listing
# get the last word of any lines that start with 'd':
 > folderLines <- items[substr(items, 1, 1)=='d']
# get the directory names and create a new data frame:
 > dirs <- unlist(lapply(strsplit(folderLines, " "), function(x){x[length(x)]}))
 > dates <- data.frame(dirname=unlist(strsplit(dirs, "\r")))
 

Next, we need to known the h/v position of the MODIS blocks. For example, we want to generate EVI images for the whole area of British Columbia. This area covers 9 MODIS tiles: h09v03, h09v04, h10v02, h10v03, h10v04, h11v02, h11v03, h12v02 and h12v03.

Figure: MODIS tiling system in MRT. In order to see which blocks you need to cover your area of interest, use the MRTWeb tool, this MODIS Sinusoidal Grid shape file, and/or this Google Earth file.

Each MODIS tile has an unique name. We can do a directory listing and get the full names of the tiles on the FTP by combining the getURL and grep methods:

 
 > getlist <- strsplit(getURL(paste(MOD13A3, dates$dirname[[1]], "/", sep=""), 
+      .opts=curlOptions(ftplistonly=TRUE\)), "\r\n")[[1]]
 > str(getlist)
chr [1:1144] "BROWSE.MOD13A3.A2000032.h00v08.005.2006271174446.1
> getlist[1101:1102]
#! [1] "MOD13A3.A2000032.h31v10.005.2006271174005.hdf"    
#! [2] "MOD13A3.A2000032.h31v10.005.2006271174005.hdf.xml"

This means that the directory "ftp://e4ftl01u.ecs.nasa.gov/MOLT/MOD13A3.005/2000.02.01/" contains a total of 1144 files. We wish to obtain only the names of the HDF files (9 tiles) for our area of interest. These can be obtained using the grep method:

 
 > BLOCK1 <- getlist[grep(getlist, pattern="MOD13A3.*.h09v03.*.hdf")[1]]
 > BLOCK2 <- getlist[grep(getlist, pattern="MOD13A3.*.h09v04.*.hdf")[1]]
...
 > BLOCK9 <- getlist[grep(getlist, pattern="MOD13A3.*.h12v02.*.hdf")[1]]
 > BLOCK1
 [1] "MOD13A3.A2000032.h09v03.005.2006271173514.hdf"
 

note that we look only for the first (HDF) file. The second file with the same name is carries the production metadata (much about the HDF tile can be read already from the file name; see the MODIS Naming Conventions).

Next, we can download each tile using the download.file method, and with help of the wget package (you need to download the wget exe to your windows system directory, otherwise you will not be able to download the tiles from R):

 
 > download.file(paste(MOD13A3a, dates$dirname[[i]], "/", BLOCK1,sep=""), 
+   destfile=paste(getwd(), "/", BLOCK1, sep=""), mode='wb', method='wget')
 --2009-03-01 18:21:37--  ftp://anonymous:*password*@e4ftl01u.ecs.nasa.gov/MOLT/...
           => `D:/PineBeetleBC/MODIS/MOD13A3.A2000032.h09v03.005.2006271173514.hdf'
 Resolving e4ftl01u.ecs.nasa.gov... 152.61.4.83
 Connecting to e4ftl01u.ecs.nasa.gov|152.61.4.83|:21... connected.
 Logging in as anonymous ... Logged in!
 ==> SYST ... done.    ==> PWD ... done.
 ==> TYPE I ... done.  ==> CWD /MOLT/MOD13A3.005/2000.02.01 ... done.
 ==> SIZE MOD13A3.A2000032.h09v03.005.2006271173514.hdf ... 3523318
 ==> PASV ... done.    ==> RETR MOD13A3.A2000032.h09v03.005.2006271173514.hdf ... done.
 Length: 3523318 (3.4M)
 

By looping these operations, you can download any tile of interest (see also the original R script used).

Mosaicking and resampling

Once we have downloaded all tiles of interest, we can mosaic them into a single image. To do this, we use the MODIS Resampling Tool, that you should have already installed on your machine. We first generate a parameter file (this is simply a list of tiles that need to be glued):

 
 > dirname1 <- sub(sub(pattern="\\.", replacement="_", dates$dirname[[i]]), 
+    pattern="\\.", replacement="_", dates$dirname[[i]])  # web-friendly filename;
# mosaic the blocks:
 > mosaicname <- file(paste(MRT, "TmpMosaic.prm", sep=""), open="wt")
 > write(paste(workd, BLOCK1, sep=""), mosaicname)
 > write(paste(workd, BLOCK2, sep=""), mosaicname, append=T)
...
 > write(paste(workd, BLOCK9, sep=""), mosaicname, append=T)
 > close(mosaicname)
 

and then run the MRT mosaicking tool using this parameter file:

 
# generate temporary mosaic:
 > shell(cmd=paste(MRT, 'mrtmosaic -i ', MRT, 'TmpMosaic.prm 
+        -s "0 1 0 0 0 0 0 0 0 0 0" -o ', workd, 'TmpMosaic.hdf', sep=""))
 

where -s "0 1 0 0 0 0 0 0 0 0 0" is the definition of the spectral selection (2nd band is the EVI image; see also the complete layer specification), and 'TmpMosaic.hdf' is the temporary mosaic HDF image.

Next, we need to generate a MRT parameter file that can be use to resample the image to the target coordinate system NAD83 / BC Albers:

 
 > filename = file(paste(MRT, "mrt", dirname1, ".prm", sep=""), open="wt")
 > write(paste('INPUT_FILENAME = ', workd, 'TmpMosaic.hdf', sep=""), filename) 
 > write('SPECTRAL_SUBSET = ( 1 )', filename, append=TRUE\)
 > write('  ', filename, append=TRUE\)
 > write('SPATIAL_SUBSET_TYPE = OUTPUT_PROJ_COORDS', filename, append=TRUE\)
 > write('  ', filename, append=TRUE\)
 > write('SPATIAL_SUBSET_UL_CORNER = ( 637278.0 1701350.0 )', filename, append=TRUE\)
 > write('SPATIAL_SUBSET_LR_CORNER = ( 1907278.0 335350.0 )', filename, append=TRUE\)
 > write('  ', filename, append=TRUE\)
 > write(paste('OUTPUT_FILENAME = ', workd, 'tmp', dirname1, '.tif', sep=""), filename, append=TRUE\)
 > write('  ', filename, append=TRUE\)
 > write('RESAMPLING_TYPE = NEAREST_NEIGHBOR', filename, append=TRUE\)
 > write('  ', filename, append=TRUE\)
 > write('OUTPUT_PROJECTION_TYPE = AEA', filename, append=TRUE\)
 > write('  ', filename, append=TRUE\)
 > write('OUTPUT_PROJECTION_PARAMETERS = ( ', filename, append=TRUE\)
 > write(' 0.0 0.0 50.0', filename, append=TRUE\)
 > write(' 58.5 -126.0 45.0', filename, append=TRUE\)
 > write(' 1000000.0 0.0 0.0', filename, append=TRUE\)
 > write(' 0.0 0.0 0.0', filename, append=TRUE\)
 > write(' 0.0 0.0 0.0 )', filename, append=TRUE\)
 > write('  ', filename, append=TRUE\)
 > write('DATUM = NAD83', filename, append=TRUE\)
 > write('  ', filename, append=TRUE\)
 > write('OUTPUT_PIXEL_SIZE = 1000', filename, append=TRUE\)
 > write('  ', filename, append=TRUE\)
 > close(filename)
 

For more info about the structure and definition of variable used in the MRT parameter files see the MRT user's guide (Land Processes DAAC, 2008).

We do not need the complete image, but only values of EVI for our area of interest. In other words, we need to resample this mosaic using the resample exe of the MRT. This will reproject and cut this mosaic using:

 
 > shell(cmd=paste(MRT, 'resample -p ', MRT, 'mrt', dirname1, '.prm', sep=""))
 

Just to check that everything is ok:

 
# Check the validity:
 > GDALinfo("tmp2000_02_01.1_km_monthly_EVI.tif")
 Closing GDAL dataset handle 0x0196f688...  destroyed ... done.
 rows        1366 
 columns     1270 
 bands       1 
 ll.x        637778 
 ll.y        1699850 
 res.x       1000 
 res.y       1000 
 oblique.x   0 
 oblique.y   0 
 driver      GTiff 
 projection  +proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45
 +lon_0=-126 +x_0=1000000 +y_0=0
 +ellps=GRS80 +datum=NAD83 +units=m
 +no_defs 
 file        tmp2000_02_01.1_km_monthly_EVI.tif 
 

The final result of resampling is shown below. Note that the final image is only 3.4 MB big, although we have downloaded almost 100 MB of data.

Figure: Final output of resampling (EVI image for july 2008), as visualized in SAGA GIS.

To analyze such time-series data (in this case 108 EVI images) you should consider using some of the free-GIS tools e.g. ILWIS GIS, GRASS GIS, or similar.

FWTools / MODIS Web services

An alternative approach to read HDF images is by using the FWTools (courtesy of by Frank Warmerdam). With FWTools you can also automate reading and resampling of MODIS bands (most probably you will first need to export the HDF bands to some standard GDAL supported format e.g. geotiff, then resample or mosaic the images). As the input coordinate system you can try using:

"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"

To work around to the Modis Reprojection Tool is to utilize FWTools to generate VRT (virtual mosaic) files by using the gdalbuildvrt functionality e.g.:

 
# create a VRT and subset with -te xmin ymin xmax ymax:
 > system(paste("gdalbuildvrt -te 755519 6618879 1606316 7809995 
+    -input_file_list ",txt.dir, "listTemp.txt ", txt.dir, "temp.vrt", sep="")) 
 > system(paste("gdalwarp -t_srs EPSG:2393 ",txt.dir, "temp.vrt ", 
+    txt.dir,"temp2.vrt", sep=""))
# reproject in the Finnish coordinates system:
 > system(paste("gdal_translate -of GTiff -projwin 3056628 7788872 
+    3810893 6608998 ",txt.dir, "temp2.vrt ", 
+    txt.dir, "MODIS_NDVI_16day_250m_", dates[i], ".tif", sep="")) 
# convert to GTiff (the projwin is not required)
 > unlink(paste(txt.dir, "temp.vrt",sep=""))
 > unlink(paste(txt.dir, "temp2.vrt",sep=""))
 > unlink(paste(txt.dir, "listTemp.txt",sep=""))
# erase the temporary files
 

This is an example suggested by Alex Villers.

There is also an (experimental) MODIS web service that provides a number of functions to fetch information about the MODIS data, e.g. getbands("MOD13Q1") will return the list of bands available for a product, and get a subset of the data. You can also download a subset of the MODIS products via a web-interface.

References


<Rating> Rate this article: 1 Poor and in error 2 Far from ready 3 Needs improvements 4 Minor corrections 5 Very useful!! </Rating>

Personal tools