Import MODIS data in GRASS using r.in.gdal

NASA offers access to its MODIS and ASTER data sets through Reverb|Echo. The data comes in HDF format and uses the Sinusoidal grid tiling system. If your gdal is compiled with HDF4 support (use ./configure –with-hdf4), you can use gdal, or any software that uses gdal, to open the downloaded MODIS tiles directly. For example in QGIS as explained here or in GRASS GIS.

In GRASS you can use the r.in.gdal function. The first step is to create a location and mapset with the Sinusoidal projection used for the MODIS data. The proj4 definition is +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs. Or you can create a new location and mapset using the MODIS Sinusoidal GIS shape file from here (this is where I got this proj4 definition from).

The next, tricky, part is that HDF4 data files can contain multiple data arrays or bands and we need to specify which of these bands we want to extract. You can use the gdalinfo function for that.

gdalinfo MOD14A1.A2012313.h22v10.hdf
...
Subdatasets:
  SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MOD14A1.A2012313.h22v10.hdf":MODIS_Grid_Daily_Fire:FireMask
  SUBDATASET_1_DESC=[8x1200x1200] FireMask MODIS_Grid_Daily_Fire (8-bit unsigned integer)
  SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:"MOD14A1.A2012313.h22v10.hdf":MODIS_Grid_Daily_Fire:QA
  SUBDATASET_2_DESC=[8x1200x1200] QA MODIS_Grid_Daily_Fire (8-bit unsigned integer)
  SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"MOD14A1.A2012313.h22v10.hdf":MODIS_Grid_Daily_Fire:MaxFRP
  SUBDATASET_3_DESC=[8x1200x1200] MaxFRP MODIS_Grid_Daily_Fire (32-bit integer)
  SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:"MOD14A1.A2012313.h22v10.hdf":MODIS_Grid_Daily_Fire:sample
  SUBDATASET_4_DESC=[8x1200x1200] sample MODIS_Grid_Daily_Fire (16-bit unsigned integer)
...

In the example above, only the relevant part of the output of gdalinfo is given; the bands included and their names. The actual information is much longer so you’ll have to look through the output to find the relevant section. To import one of these bands in GRASS, you need the name as specified by gdalinfo.

For example, you have the file MOD14A1.A2012313.h22v10.hdf and you want to extract the FireMask band. That means you need to name of the first band, which is “HDF4_EOS:EOS_GRID:”MOD14A1.A2012313.h22v10.hdf”:MODIS_Grid_Daily_Fire:FireMask”

You can do this manually, but that can become a bit tedious. Easier is to use the command line utilities grep and sed, like in the example below:

inputmap=`gdalinfo MOD14A1.A2003097.h21v07.hdf | grep 'SUBDATASET_1_NAME=' | sed -e 's/SUBDATASET_1_NAME=//'`
outputmap=`echo $inputmap |  awk -F\" '{ print $2 }' | awk -F. '{ print $1"_"$2"_"$3 }' | sed -e 's/A//'`
r.in.gdal $inputmap out=$outputmap

The first line extracts the name of SUBDATASET_1_NAME. The second creates the name of the imported layer in your GRASS database (MOD141_A2003097_h21v07 in the example above). The third line is the actual r.in.gdal command, which imports the first band of the file in the current mapset.

a little bit more involving. However, the advantage is that it is relatively easy to create a script based on these few lines for batch importing data.

Check also out the GRASS GIS Wiki, where you can find a few alternative ways to import the MODIS data layers. You might also be interested to read this article on spatial-analyst.net. And of course, NASA offers the Modis Reprojection Tool (MRT), to reproject and mosaic the MODIS hdf data tiles.

Advertisements

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