Thu 27 June 2013

Filed under Blog

Tags MODIS GDAL tips

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:

  1. Downloading the data
  2. Stitching together the data
  3. Reprojecting
  4. 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
Comment

Sun 23 June 2013

Filed under Blog

Tags GDAL tips

Let's say one wants to find out a box around a given latitude/longitude point. How to do this without GIS software using an artisan UNIX approach? Hell, if it works for coffee! The steps are as follows:

  1. Convert the centre point into UTM so we have meters rather ...
Read More

Thu 23 May 2013

Filed under Blog

Tags python gdal tips

This note details how to use parallel processing in python accessing shared memory. The usage case is when you have a process that reads in a lot of data, and where the data can be processed in parallel, such as when reading and processing images/stacks of images in a ...

Read More

Wed 08 May 2013

Filed under Blog

Tags unix tips meteo

A common requirement when dealing with processing large datasets over multiple networked machines is to have a local staging space: copy the data on the local disk to improve access speed and not to bog down a NFS server when all different processes start accessing different files all at once ...

Read More

Tue 11 December 2012

Filed under Blog

Tags python tips meteo

Many networks of automated meteo stations are around. In some places, the data are easy to access and therefore use. In some others, you have to go through a number of hoops. In Andalusia, it's a mixed bag: a number of independent networks are available for agricultural and environmental ...

Read More

EO & DA Ramblings © J Gomez-Dans Powered by Pelican and Twitter Bootstrap. Icons by Font Awesome and Font Awesome More