Stitching MODIS burned area datasets

Thu 27 June 2013

Filed under Blog

MODIS products are usually provided as data granules, representing the magnitude of interest (plus a number of different extra layers of metadata, quality assurance flags, etc) over a given temporal period for an area typically extending 1200 by 1200 km in the MODIS sinusoidal projection, around 10x10 degrees. Typically, your study site might extend over several tiles, or just happen to lie between two tiles. Consider Portugal, split (mostly, there are islands and other stuff, but we typically ignore them ;D) into two tiles, h17v04 and h17v05. In this post, we look at:

Downloading the data
Stitching together the data
Reprojecting
Applying a mask

Downloading the data

We are going to use the get_modis script (you can get it from here). It is a simple Python script that automatically downloads MODIS data. We shall download data for the year 2004. What we need is the following command

./get_modis.py --verbose --platform=MOTA --product=MCD45A1.051 --tile=h17v04 --year=2004 --output=.

where

platform=MOTA: This is the "platform". In this case, the MCD45A1 product uses data from TERRA and AQUA, so its code is MOTA.
product=MCD45A1.051: We want the MCD45A1 product, collection 5.1 (just released!).
tile=h17v04: A tile for the Northern Iberian Peninsula.
year=2004: The year of interest
output=.: Use the current directory to save the downloaded files.

Because we provided the --verbose flag, we get a running commentary on how all this proceeds. For the other tile (h17v05), this is the output

   2013-06-27 13:31:58,898 Getting MCD45A1.A2004001.h17v05.051.2013070231310.hdf.....
   2013-06-27 13:32:01,840 Done!
   2013-06-27 13:32:04,604 Getting MCD45A1.A2004032.h17v05.051.2013062125251.hdf.....
   2013-06-27 13:32:06,706 Done!
   2013-06-27 13:32:08,207 Getting MCD45A1.A2004061.h17v05.051.2013071001320.hdf.....
   2013-06-27 13:32:10,474 Done!
   2013-06-27 13:32:12,046 Getting MCD45A1.A2004092.h17v05.051.2013062142742.hdf.....
   2013-06-27 13:32:14,177 Done!
   2013-06-27 13:32:15,909 Getting MCD45A1.A2004122.h17v05.051.2013062153940.hdf.....
   2013-06-27 13:32:17,987 Done!
   2013-06-27 13:32:19,521 Getting MCD45A1.A2004153.h17v05.051.2013062164329.hdf.....
   2013-06-27 13:32:21,770 Done!
   2013-06-27 13:32:24,424 Getting MCD45A1.A2004183.h17v05.051.2013062175918.hdf.....
   2013-06-27 13:32:26,220 Done!
   2013-06-27 13:32:28,791 Getting MCD45A1.A2004214.h17v05.051.2013062184842.hdf.....
   2013-06-27 13:32:30,599 Done!
   2013-06-27 13:32:32,367 Getting MCD45A1.A2004245.h17v05.051.2013062201148.hdf.....
   2013-06-27 13:32:34,460 Done!
   2013-06-27 13:32:36,893 Getting MCD45A1.A2004275.h17v05.051.2013062212709.hdf.....
   2013-06-27 13:32:39,459 Done!
   2013-06-27 13:32:41,015 Getting MCD45A1.A2004306.h17v05.051.2013062214733.hdf.....
   2013-06-27 13:32:43,370 Done!
   2013-06-27 13:32:45,977 Getting MCD45A1.A2004336.h17v05.051.2013062215105.hdf.....
   2013-06-27 13:32:50,383 Done!
   2013-06-27 13:32:50,383 Completely finished downlading all there was


Stitching the data
-------------------

MCD45A1 granules do have a number of different datasets. You can have a peek inside using ``gdalinfo``:

::

    gdalinfo MCD45A1.A2004214.h17v04.051.2013062164342.hdf
    [ Lots and lots and lots of metadata]
    Subdatasets:
     SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MCD45A1.A2004214.h17v04.051.2013062164342.hdf":MOD_GRID_Monthly_500km_BA:burndate
     SUBDATASET_1_DESC=[2400x2400] burndate MOD_GRID_Monthly_500km_BA (16-bit integer)
     SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:"MCD45A1.A2004214.h17v04.051.2013062164342.hdf":MOD_GRID_Monthly_500km_BA:ba_qa
     SUBDATASET_2_DESC=[2400x2400] ba_qa MOD_GRID_Monthly_500km_BA (8-bit unsigned integer)
     SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"MCD45A1.A2004214.h17v04.051.2013062164342.hdf":MOD_GRID_Monthly_500km_BA:npasses
     SUBDATASET_3_DESC=[2400x2400] npasses MOD_GRID_Monthly_500km_BA (8-bit unsigned integer)
     SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:"MCD45A1.A2004214.h17v04.051.2013062164342.hdf":MOD_GRID_Monthly_500km_BA:nused
     SUBDATASET_4_DESC=[2400x2400] nused MOD_GRID_Monthly_500km_BA (8-bit unsigned integer)
     SUBDATASET_5_NAME=HDF4_EOS:EOS_GRID:"MCD45A1.A2004214.h17v04.051.2013062164342.hdf":MOD_GRID_Monthly_500km_BA:direction
     SUBDATASET_5_DESC=[2400x2400] direction MOD_GRID_Monthly_500km_BA (8-bit unsigned integer)
     SUBDATASET_6_NAME=HDF4_EOS:EOS_GRID:"MCD45A1.A2004214.h17v04.051.2013062164342.hdf":MOD_GRID_Monthly_500km_BA:surfacetype
     SUBDATASET_6_DESC=[2400x2400] surfacetype MOD_GRID_Monthly_500km_BA (8-bit unsigned integer)
     SUBDATASET_7_NAME=HDF4_EOS:EOS_GRID:"MCD45A1.A2004214.h17v04.051.2013062164342.hdf":MOD_GRID_Monthly_500km_BA:gap_range1
     SUBDATASET_7_DESC=[2400x2400] gap_range1 MOD_GRID_Monthly_500km_BA (16-bit unsigned integer)
     SUBDATASET_8_NAME=HDF4_EOS:EOS_GRID:"MCD45A1.A2004214.h17v04.051.2013062164342.hdf":MOD_GRID_Monthly_500km_BA:gap_range2
     SUBDATASET_8_DESC=[2400x2400] gap_range2 MOD_GRID_Monthly_500km_BA (16-bit unsigned integer)
   Corner Coordinates:
   Upper Left  (    0.0,    0.0)
   Lower Left  (    0.0,  512.0)
   Upper Right (  512.0,    0.0)
   Lower Right (  512.0,  512.0)
   Center      (  256.0,  256.0)

We are intersted (initially) in the first dataset burndate. We want to, for each month, stich together the two tiles, or mosaic them. To make it easy and efficient, we will do it in two stages using virtual datasets. These are just text files that refer to the original datasets and allow fast modifications. Fist, we create VRTs for the burndate layers for the two tiles. In BASH, this can be done easily with a loop

for fich in *h17v0*.hdf
do
    gdal_translate -of VRT 'HDF4_EOS:EOS_GRID:"'${fich}'":MOD_GRID_Monthly_500km_BA:burndate' \
        `basename ${fich} .hdf`.vrt
done

We now have a number of files with the extension .vrt. You can read their contents, these are just text files (XML actually) and point to the original datasets. We now use gdalwarp to mosaic them, and to just get the Portugal mask. To this end, we use the world.shp of APRS world. Since this shapefile does not have a projection, we set it to WGS84 Lat/Long:

ogr2ogr -a_srs "EPSG:4326" world_wgs84.shp world.shp

This gives a new file, world_wgs84.shp with a projection in it. We then just use gdalwarp to reproject and cut, but for completenes, we shall also produce data for Spain.

for fich in MCD45A1.A2004*h17v04*vrt
do
    other=`echo $fich | cut -d "." -f 1,2,`.h17v05.*vrt
    output_portugal=`echo $fich | cut -d "." -f 1,2`.Portugal.tif
    output_spain=`echo $fich | cut -d "." -f 1,2`.Spain.tif
    gdalwarp -of GTiff -co "COMPRESS=LZW" \
        -cutline world_wgs84.shp -cwhere "NAME='PORTUGAL'" \
        ${fich} ${other} ${output_portugal}
    gdalwarp -of GTiff -co "COMPRESS=LZW" \
        -cutline world_wgs84.shp -cwhere "NAME='SPAIN'" \
        ${fich} ${other} ${output_spain}
done

We now have a number of GeoTIFF files, but these are in the MODIS projection. Let's say that we wanted them in UTM (UTM30N/ED50 for Portugal and UTM29N/ED50 for Spain). We'd go about it using the relevant EPSG codes as:

for fich in MCD45A1.A2004*h17v04*vrt
do
    other=`echo $fich | cut -d "." -f 1,2,`.h17v05.*vrt
    output_portugal=`echo $fich | cut -d "." -f 1,2`.Portugal.UTM.tif
    output_spain=`echo $fich | cut -d "." -f 1,2`.Spain.UTM.tif
    gdalwarp -of GTiff -co "COMPRESS=LZW" -t_srs "EPSG:23030" \
        -cutline world_wgs84.shp -cwhere "NAME='PORTUGAL'" \
        ${fich} ${other} ${output_portugal}
    gdalwarp -of GTiff -co "COMPRESS=LZW" -t_srs "EPSG:23029" \
        -cutline world_wgs84.shp -cwhere "NAME='SPAIN'" \
        ${fich} ${other} ${output_spain}
done

Stitching MODIS burned area datasets

Downloading the data

Comments