PlanetFlux

Processing MODIS Swath data with the HEG tool

2012-10-15T11:32:00.001-04:00

I recently found a bug/weakness in the NASA's HEG tool which can be used to convert MODIS swath data (such as the MOD06_L2 cloud product) into a regular gridded file. The library requires several environment parameters to be set correctly before running successfully. For example:

MRTDATADIR=.../heg/data
PGSHOME=.../heg/TOOLKIT_MTD

When those parameters are set incorrectly, the program does not give any indication of an error and produces an output file full of gibberish. For example, here is an image of a sample MOD06 scene (MOD06_L2.A2008062.1545.051.2010324234929.hdf) of "cloud particle effective radius" from the swath HDF file acquired from NASA (visualized using HDFview):

The black areas are not cloudy and the colorful areas represent cloudy areas with different properties. It intuitively looks like a reasonable image of cloudiness. After conversion to a gridded HDF file using the HEG program swtif and incorrectly set environment variables, this is what the image looks like:

Obviously there are some conversion issues (this image looks completely different), but after running the program there was no indication that an error had occurred (either in STOUT or in the swtif.log file). Only after a few hours of trial and error (trying different projections, etc.) did I realize what was causing the problem. While I realize accurately following the instructions is the responsibility of the user, it is still a bit dangerous for a program to output erroneous results when important libraries cannot be found. It would be trivial to add a quick check that the libraries are readable before writing the output file.

I've also sent this report to the HEG folks, but I'm posting this here in case others find themselves in a similar situation.

Compiling GRASS on a cluster running SGI ProPack for Linux

2012-10-05T07:39:00.001-04:00

Some notes on installing GRASS (and a few dependencies) and NCO from source in a personal directory (no su privileges) on a cluster running SGI ProPack for Linux . I mostly followed the instructions here.

Set install path and add it to PATH

set IPATH=(INSTALL LOCATION)
setenv PATH $IPATH\bin:$PATH
echo 'PATH=$IPATH/bin:$PATH' >> ~/.profile

Load necessary modules

module load gcc hdf5 netcdf/4.1.3/gcc/mpt mpi tcl-tk/8.5.11 udunits/2.1.19 szip/2.1/gcc

Compile proj4

svn checkout http://svn.osgeo.org/metacrs/proj/trunk/proj
cd proj
./configure --prefix=$IPATH
make
make install

Compile GEOS

wget http://download.osgeo.org/geos/geos-3.3.5.tar.bz2
tar xvf geos-3.3.5
cd geos-3.3.5
./configure --prefix=$IPATH
make -j8
make install

Compile HDF4 (for some reason the built in version wasn't working)

wget http://www.hdfgroup.org/ftp/HDF/HDF_Current/src/hdf-4.2.8.tar.gz
tar -xvf hdf-4.2.8.tar.gz
cd hdf-4.2.8
./configure --prefix=$IPATH --with-szlib=/nasa/szip/2.1 --disable-netcdf --disable-fortran --enable-shared

Compile GDAL

wget http://download.osgeo.org/gdal/gdal-1.9.1.tar.gz
tar xvf gdal-1.9.1
cd gdal-1.9.1
./configure FC=mpif90 CC=mpicc CXX=mpicxx --prefix=$IPATH --with-netcdf=/nasa/sles11/netcdf/4.1.3/gcc/mpt/bin --with-geos=$IPATH/bin/geos-config --with-static-proj4=$IPATH/ --with-hdf5=/nasa/hdf5/1.8.3/gcc/mpt/lib  --with-hdf4=/nobackupp1/awilso10/software/lib make -j8
make install

Checkout the current version of GRASS via SVN:

svn checkout https://svn.osgeo.org/grass/grass/branches/releasebranch_6_4 grass64_release
./configure --with-cxx --with-proj-share=$IPATH/share/proj --with-proj-libs=$IPATH/lib --with-proj-includes=$IPATH/include --without-fftw --with-geos=$IPATH/lib --with-blas-includes=/nasa/intel/mkl/10.0.011/include/ --with-blas-libs=/nasa/intel/mkl/10.0.011/lib --prefix=$IPATH
make
make install

Also install CDO tools to process netcdf files

wget https://code.zmaw.de/attachments/download/3873/cdo-current.tar.gz
tar -xvf cdo-current.tar.gz 
./configure --prefix=$IPATH --with-netcdf=/nasa/sles11/netcdf/4.1.3/gcc/mpt --with-hdf5=/nasa/hdf5/1.8.3/gcc/mpt
make
make install

Compile Rmpi
This took a lot of time to figure out because most versions of MPI available on this system do not seem to work with Rmpi...

wget http://cran.r-project.org/src/contrib/Rmpi_0.6-1.tar.gz
module load mpi-intel
export MPIHOME=/nasa/intel/impi/4.0.0.028/
CC=mpicc R CMD INSTALL --configure-args="--with-Rmpi-type=OPENMPI  --with-Rmpi-libpath=$MPIHOME/lib64 --with-Rmpi-include=$MPIHOME/include64" ~/Rmpi_0.6-1.tar.gz

Photographing the Miniature Forests of Cape Horn

2012-09-07T09:42:00.000-04:00

In January 2011 I spent five weeks in the Chilean Patagonia working to photograph the incredible bryophyte diversity in the "Miniature Forests of Cape Horn."

Bernard Goffinet on the slopes of Mt. Robelo

The first of two books from this project is now available from the University of North Texas Press (and Amazon).

Miniature Forests of Cape Horn: Ecotourism with a Hand Lens

Bernard Goffinet, Ricardo Rozzi, Lily Lewis, William Buck, Francisca Massardo

A video documenting the ongoing socio-ecological research there is also available:

Regreso a la Madriguera from Cerebro Producciones on Vimeo. Here's the description in English
"All the beauty of the National Park in Tierra del Fuego, Omora serves as the setting for this story documentary which exposes a new way to teach ecology, from an ethical and humane perspective. Philosophers and biologists at the University of Texas and the Magellan lead us to the wealth of one of the world's southernmost forests."

What your projection says about you

2012-06-10T10:57:00.000-04:00

Great cartoon from XKCD.com

Speeding up R with Intel's Math Kernel Library (MKL)

2012-05-02T10:24:00.001-04:00

I did some comparisons of the generic BLAS with Intel's MKL (both sequential and parallel) on a Dell PowerEdge 610 server with dual hyperthreading 6-core 3.06GHz Xeon X5675 processors. Here are the results from an R benchmarking script (Normal R indicates the generic BLAS, sMKL is the sequential (single core Intel MKL, and pMKL is the parallel Intel MKL using all 24 threads available on this system). Times are in seconds, lower is better.

R Benchmark 2.5

===============

Number of times each test is run__________________________: 3

I. Matrix calculation

--------------------- Normal R sMKL pMKL

Creation, transp., deformation of a 2500x2500 matrix (sec): 0.592 0.583 0.585

2400x2400 normal distributed random matrix ^1000____ (sec): 0.425 0.411 0.427

Sorting of 7,000,000 random values__________________ (sec): 0.787 0.778 0.777

2800x2800 cross-product matrix (b = a' * a)_________ (sec): 11.543 1.875 0.283

Linear regr. over a 3000x3000 matrix (c = a \ b')___ (sec): 5.367 0.910 0.214

Trimmed geom. mean (2 extremes eliminated): 1.358 0.743 0.414

II. Matrix functions

--------------------

FFT over 2,400,000 random values____________________ (sec): 0.422 0.451 0.435

Eigenvalues of a 640x640 random matrix______________ (sec): 0.949 0.443 0.414

Determinant of a 2500x2500 random matrix____________ (sec): 4.864 0.967 0.352

Cholesky decomposition of a 3000x3000 matrix________ (sec): 4.131 0.865 0.179

Inverse of a 1600x1600 random matrix________________ (sec): 4.011 0.751 0.277

Trimmed geom. mean (2 extremes eliminated): 2.505 0.667 0.343

III. Programmation

------------------

3,500,000 Fibonacci numbers calculation (vector calc)(sec): 0.787 0.824 0.841

Creation of a 3000x3000 Hilbert matrix (matrix calc) (sec): 0.456 0.465 0.431

Grand common divisors of 400,000 pairs (recursion)__ (sec): 2.196 2.386 1.927

Creation of a 500x500 Toeplitz matrix (loops)_______ (sec): 0.616 0.612 0.596

Escoufier's method on a 45x45 matrix (mixed)________ (sec): 0.470 0.425 0.447

Trimmed geom. mean (2 extremes eliminated): 0.611 0.617 0.607

--------------------------------------------

Total time for all 15 tests_________________________ (sec): 37.62 12.76 8.18

Overall mean (sum of I, II and III trimmed means/3)_ (sec): 1.28 0.67 0.44

So you can see there are some significant gains, especially for the slowest tasks (see bolded items above). For example, the parallel MKL resulted in a 40x speedup of the cross-product calculation. Across all jobs the parallel MKL version was ~4.5x faster than the generic BLAS on this system. When running models that take days to fit, that is significant!

Here's how I installed it on Ubuntu 12.04:

Download and install the Intel MKL from here.
First install some libraries needed by R:sudo apt-get install libreadline6 libreadline6-dev xserver-xorg xserver-xorg-dev gfortran-devel libxt-dev
Download the latest R, I used v2.15, available here
Compile R from source using something like this (detailed instructions here):
export MKL_LIB_PATH=/opt/intel/composer_xe_2011_sp1.7.256/mkl/lib/intel64export optim_flags="-O3 -funroll-loops -march=native"# from http://cran.r-project.org/doc/manuals/R-admin.html#MKL=" -L${MKL_LIB_PATH} -Wl,--start-group ${MKL_LIB_PATH}/libmkl_gf_lp64.a \ ${MKL_LIB_PATH}/libmkl_gnu_thread.a ${MKL_LIB_PATH}/libmkl_core.a \ -Wl,--end-group -lgomp -lpthread"sudo ./configure --enable-R-shlib --with-blas="$MKL" CC="gcc" CFLAGS="$optim_flags" CXX="g++" CXXFLAGS="$optim_flags" F77="gfortran" FFLAGS="$optim_flags" FC="gfortran" FCFLAGS="$optim_flags"

Installing R's maps package on Ubuntu

2012-04-16T14:21:00.002-04:00

I recently ran into trouble trying to install the R maps package on Ubuntu 10.04. Here's the error I was getting:

** arch -
gcc -std=gnu99 -O3 -pipe -g Gmake.c -o Gmake
Gmake.c: In function ‘get_lh’:
Gmake.c:111: warning: cast from pointer to integer of different size
Gmake.c:113: warning: cast from pointer to integer of different size
Gmake.c: In function ‘main’:
Gmake.c:211: warning: cast from pointer to integer of different size
Gmake.c:214: warning: cast from pointer to integer of different size
Gmake.c:217: warning: cast from pointer to integer of different size
Gmake.c:219: warning: cast from pointer to integer of different size
Gmake.c:221: warning: cast from pointer to integer of different size
Gmake.c:224: warning: cast from pointer to integer of different size
Gmake.c:227: warning: cast from pointer to integer of different size
gcc -std=gnu99 -O3 -pipe -g Lmake.c -o Lmake
Lmake.c: In function ‘main’:
Lmake.c:223: warning: cast from pointer to integer of different size
Lmake.c:228: warning: cast from pointer to integer of different size
Lmake.c:230: warning: cast from pointer to integer of different size
Lmake.c:232: warning: cast from pointer to integer of different size
Lmake.c:235: warning: cast from pointer to integer of different size
Converting world to world2
f convert.awk < world.line > world2.line
/bin/bash: f: command not found
make: [world2.line] Error 127 (ignored)
make county.L state.L usa.L nz.L world.L world2.L italy.L france.L state.vbm.L state.carto.L
make[1]: Entering directory `/tmp/Rtmp4lXbpt/R.INSTALL375371de/maps/src'
./Lmake 0 s b county.line county.linestats ../inst/mapdata/county.L
./Lmake 0 s b state.line state.linestats ../inst/mapdata/state.L
./Lmake 0 s b usa.line usa.linestats ../inst/mapdata/usa.L
./Lmake 0 s b nz.line nz.linestats ../inst/mapdata/nz.L
./Lmake 0 s b world.line world.linestats ../inst/mapdata/world.L
./Lmake 0 s b world2.line world2.linestats ../inst/mapdata/world2.L
Cannot read left and right at line 1
make[1]: *** [world2.L] Error 1
make[1]: Leaving directory `/tmp/Rtmp4lXbpt/R.INSTALL375371de/maps/src'
make: *** [ldata] Error 2
ERROR: compilation failed for package ‘maps’
* removing ‘/usr/local/lib/R/site-library/maps’
ERROR: dependency ‘maps’ is not available for package ‘ramps’
* removing ‘/usr/local/lib/R/site-library/ramps’

The downloaded packages are in
‘/tmp/RtmpsAvFte/downloaded_packages’

After banging my head on it for awhile, I eventually found that if you

download the tar.gz from CRAN
change line 55 of the Makefile to:awk -f convert.awk < world.line > world2.line

it will install fine.

Slow R graphics via SSH

2012-03-30T11:00:00.004-04:00

I do most of my R work on a server via SSH. Viewing graphics (i.e. plot()) can be excruciatingly slow through the SSH pipe, but it's easy to make it fast. Simply run:
X11.options(type="Xlib")
once in the R session.

Installing grass extensions on Ubuntu

2012-02-14T12:23:00.001-05:00

I recently struggled to install extensions using the version of GRASS installed from the UbuntuGIS repositories. For example, I was trying to install i.landsat.toar to convert landsat's DNs to top of atmosphere radiance values using the following command:

 g.extension -s --verbose extension=i.landsat.toar svnurl=https://svn.osgeo.org/grass/grass-addons/grass6

That gave the following error during compilation:

 Using WGET for downloading data.  
 Fetching <i.landsat.toar> from GRASS-Addons SVN (be patient)...  
 A  i.landsat.toar/landsat_set.c  
 A  i.landsat.toar/local_proto.h  
 A  i.landsat.toar/main.c  
 A  i.landsat.toar/description.html  
 A  i.landsat.toar/landsat.c  
 A  i.landsat.toar/earth_sun.c  
 A  i.landsat.toar/landsat.h  
 A  i.landsat.toar/landsat_met.c  
 A  i.landsat.toar/Makefile  
 A  i.landsat.toar/earth_sun.h  
  U  i.landsat.toar  
 Checked out revision 50813.  
 Compiling <i.landsat.toar>...  
 test -d OBJ.x86_64-pc-linux-gnu || mkdir -p OBJ.x86_64-pc-linux-gnu  
 gcc -I/build/buildd/grass-6.4.1/dist.x86_64-pc-linux-gnu/include -Wall -g -O    -DPACKAGE=\""grassmods"\" -I/build/buildd/grass-6.4.1/dist.x86_64-pc-linux-gnu/include -o OBJ.x86_64-pc-linux-gnu/earth_sun.o -c earth_sun.c  
 gcc -I/build/buildd/grass-6.4.1/dist.x86_64-pc-linux-gnu/include -Wall -g -O    -DPACKAGE=\""grassmods"\" -I/build/buildd/grass-6.4.1/dist.x86_64-pc-linux-gnu/include -o OBJ.x86_64-pc-linux-gnu/landsat.o -c landsat.c  
 landsat.c:4:23: fatal error: grass/gis.h: No such file or directory  
 compilation terminated.  
 make: *** [OBJ.x86_64-pc-linux-gnu/landsat.o] Error 1  
 ERROR: Compilation failed, sorry. Please check above error messages.

After exploring the information here, here, and especially here, I finally discovered a simple solution that seems to work. I don't think this has been reported elsewhere.

First download the source code using svn as listed on the add-ons page. For example:

 svn co https://svn.osgeo.org/grass/grass-addons/grass6/imagery/i.landsat.toar

then run:

 sudo make MODULE_TOPDIR=/usr/lib/grass64 GRASS_HOME=/usr/lib/grass64 ARCH_DISTDIR=/usr/lib/grass64

It does return a few errors, but the addons are installed and work correctly (as far as I can tell).

UPDATE (2012-11-9): Following the release of GRASS 6.4.2, i.landsat.acca and i.landsat.toar are supposed to be included in core GRASS. However, I just upgraded to 6.4.2 via the UbuntuGIS ppa and neither program seem to be included. Strangely, i.landsat.toar is still available via g.extension, but i.landsat.acca is no longer in the addons SVN. By running:

 svn list -R https://svn.osgeo.org/grass/ | grep "i.landsat.acca"

I was able to find it here: grass/branches/develbranch_6/imagery/i.landsat.acca/. So just replace the download command above with this:

 svn co https://svn.osgeo.org/grass/grass/branches/develbranch_6/imagery/i.landsat.acca/

followed by the make command above and you should be back in business.

Biodiversity of Earth's Richest Plant Kingdom Under Fire

2012-01-31T09:37:00.000-05:00

The National Science Foundation just published a story on a research project I'm involved in to understand the dimensions of biodiversity. A number of my photographs were used in the story, including this time lapse movie of several scientists (including me) conducting a vegetation survey in South Africa to track changes in biodiversity.

Copy all the files in a directory to a new directory using R

2011-10-27T15:48:00.000-04:00

Someone asked me how to move a directory full of files from one place to another using R. The easiest way I've found is as follows (where "oldpath" is the existing directory and "newpath" is the new directory):

file.copy(list.files(oldpath),newpath)

Lightworks is going open source!

2011-07-01T09:19:00.000-04:00

I just became aware of Lightworks video editing software, which has been used to create several major motion pictures (including The King’s Speech).

They are taking the big (and admirable) step of making their software free and open source in order to encourage use and development (they must be planning to switch to selling services, classes, etc.). It's functionality is very similar to final cut pro (a steep learning curve, but very powerful) and it's price (free) is quite attractive. They only have a windows version available now, but there should be versions for mac and linux this year (as soon as they release the source code, people will be jumping all over this)...

Extracting reflectance data from SpectraSuite JCAMP files in R

2011-05-06T10:02:00.002-04:00

I've been using an OceanOptics USB4000 spectrometer in research on biomass accumulation and climate in South Africa.

This post presents a R function that will read in JCAMP files from the SpectraSuite software and output a vector of reflectance values. SpectraSuite can convert binary output files to JCAMP files if they weren't saved as JCAMP.

You can use this to process a folder of spectra and return a data frame with a column for each file as follows:

And here's another function that will calculate NDVI similar to several satellite products:

You can apply this to the output of spectracompile() like this:

Here's a few random spectra from various fynbos species:

Reconnecting the Lower Rio Grande

2011-04-15T13:24:00.004-04:00

The 2011 North American Nature Photographers Association College Scholarship winners recently finished their project documenting efforts to 'reconnect the Lower Rio Grande valley in Southeastern Texas. In a mere 5 days (and a few sleepless nights), they developed, planned, shot, and edited the project. A great accomplishment. Let me know what you think...

You can read their introduction on their Vimeo site, or just watch it here:

Reconnecting the Rio Grande Valley from NANPA College Students on Vimeo.

o para que usted español hablantes:

NANPA PROJECT spanish version2 from NANPA College Students on Vimeo.

One of the participants also put together an excellent documentary on the development of the project. You can see that here:

From Passion to Profession from Aaron Schmidt on Vimeo.

Big Bang, Big Boom

2011-02-02T15:38:00.005-05:00

I just ran across this incredible stop-motion video covering the last 14 billion years. Those animated murals are fantastic!

BIG BANG BIG BOOM:
an unscientific point of view on the beginning and evolution of life ... and how it could probably end.

direction and animation by BLU blublu.org
production and distribution by ARTSH.it artsh.it
soundtrack by ANDREA MARTIGNONI

Carbon counter

2010-11-12T10:37:00.011-05:00

This morning I met briefly with Bruce M. Kahn, Ph.D. who works for Deutsche Bank Climate Change Advisors. He showed me their "Carbon Counter" that presents a real-time Semi-Empirical Representation of the Temporal Variation of Total Greenhouse Gas Levels Expressed as Equivalent Levels of Carbon Dioxide (you can read the methodology here). It's on display in New York City as a huge billboard. Here's the current total:

www.know-the-number.com

Our Climate is Changing!

Please download Flash Player.

Deutsche Bank has published several reports on climate change and risk in financial markets. Here's a few lines from their website:

Three years ago, Deutsche Asset Management identified climate change as one of the mega-trends that would drive the global asset management business for the next generation and beyond. We saw that the rapidly growing level of carbon in the atmosphere meant that the world had to take action now, and that this would require massive capital investment over several decades. That in turn would produce exciting new investment opportunities from which our clients could benefit.

It makes sense that companies that make their business dealing with risk (such as banks and insurance providers) have been among the quickest to move from denying that climate change is real to accepting and planning for it.

Regarding their counter, when I saw it I wondered if the number would go down in the northern hemisphere spring (when all plants in the north start drawing down CO2), but they've removed those natural signals. Here's what the actual CO2 concentration over the last few decades looks like:

So the overall trend is definitely increasing, but there is an annual cycle that increases and decreases as the northern forests "breath" in and out every year. The counter is displaying something close to the red line on this plot, which is reasonable approach to convey that the overall trend is upward. Otherwise the people of New York might think we solved the greenhouse gas emission problem next May...

In fact, according to the Numerical Terradynamic Simulation Group who have analysed satellite data to estimate net primary productivity over the last decade, there seems to have been a decrease in plant productivity. Check out this short summary video:

Potential Solar Radiation on the Cape Of Good Hope, South Africa

2010-10-08T09:58:00.003-04:00

Incoming solar radiation drives the earth's climate system and is the ultimate source of energy on which most life on earth depends. Here's an animation of potential solar radiation on the Cape of Good Hope, South Africa:

The units are watts per square meter per day. A small bulb like those used on a Christmas tree are about one watt. So the minimum value of 375 would be nearly equivalent to having 375 of those little bulbs burning on a square meter of land for 24 hours. Now imagine the maximum value of 10,000.

These data were generated using the GRASS GIS r.sun routine using a 30m elevation model from the ASTER project.

Absolute Certainty

2010-10-07T15:43:00.002-04:00

The trouble with waiting until Absolute Certainty
(which science will never provide):

This year’s Union of Concerned Scientists' Science Idol is Raul de la Nuez of Miami, Florida. His calendar will be on the cover of the 2011 Scientific Integrity Calendar. Learn more about the calendar here.

ResearcherID

2010-09-21T15:09:00.000-04:00

I just found the ResearcherID service by Thomson Reuters that allows researchers to create a unique ID and keep track of publications (as well as automatically calculate their H-index). The Labs feature also offers some summary data of who your co-authors are (if you can't remember) and who's citing you. Check out the map of people that have cited my papers here. You have control over which publications are in your list so you don't have to rely on an automatic search. They offer a "badge" that you can embed in any html page to give readers access to an up-to-date list of your publications and citations. Here's mine:

Nature Photography Scholarship

2010-06-18T16:06:00.004-04:00

Remember the last time you talked to your family and non-scientist friends about your research and after 45 seconds their eyes glossed over. Then you pulled out some scatter plots you find really compelling, and slapped some more wet newspaper on their interest.

Now imagine instead that you pull up a photographic slide show that shows them, rather than telling them, what you've been up to. A series of well-composed photos sets the stage, leads them through your process, and draws them in. After you've piqued their interest, you might even be able to sneak in a few graphs that hammer the points home.

The North American Nature Photographer's Association is offering student scholarships to attend next year's nature photography meeting in McAllen, Texas March 7 - 13, 2011 (you don't have to be a US citizen). I got one of these last year and it was awesome. They don't pay for travel, but they cover everything else (accommodation, food, registration, etc., which are worth about $1,500). The application deadline is Novermber 8, 2010. Get more information at the NANPA website.

Download the flyer here

Shrinking R's PDF output

2010-06-17T11:28:00.000-04:00

R is great for graphics, but I've found that the PDF's R produces when drawing large plots can be extremely large. This is especially common when using spplot() to plot a large raster. I've made a 15 page PDF full of rasters that was hundreds of MB in size. Obviously I don't need all the detail (every pixel of the raster) represented in the pdf and would rather have it reduced in size somehow. So I wrote an R function to automate the following:

take an existing pdf and run ps2pdf on it as an intial compression step. Often this step is all that's needed.
split it into separate files using pdftk
Check to see if each separate page is larger than some threshold you specify (I set 5MB as the default)
If any one page is larger, rasterize the whole thing to a PNG file using ghostscript. I used the multicore package to parallelize this step, but this isn't necessary and that call could be replaced by lapply() to run them sequentially.
Put the separate pages (perhaps a mix of the original and the compressed rasters) back together.

Here's the function:

 shrinkpdf<-function(pdf,maxsize=5,suffix="_small",verbose=T){  
  require(multicore)  
   wd=getwd()  
   td=paste(tempdir(),"/pdf",sep="")  
   if(!file.exists(td)) dir.create(td)  
   if(verbose) print("Performing initial compression")  
   system(paste("ps2pdf ",pdf," ",td,"/test.pdf",sep=""))  
   setwd(td)  
   system(paste("pdftk ",td,"/test.pdf burst",sep=""))  
   files=list.files(pattern="pg_")  
   sizes=sapply(files,function(x) file.info(x)$size)*0.000001 #get sizes of individual pages  
   toobig=sizes>=maxsize  
   if(verbose)  print(paste("Resizing ",sum(toobig)," pages:  (",paste(files[toobig],collapse=","),")",sep=""))  
   mclapply(files[toobig],function(i){  
    system(paste("gs -dBATCH -dTextAlphaBits=4 -dNOPAUSE -r300 -q -sDEVICE=png16m -sOutputFile=",i,".png ",i,sep=""))  
    system(paste("convert -quality 100 -density 300 ",i,".png ",strsplit(i,".",fixed=T)[[1]][1],".pdf ",sep=""))  
    if(verbose) print(paste("Finished page ",i))  
    return()  
   })  
   if(verbose) print("Compiling the final pdf")  
   file.remove("test.pdf")  
   file.remove(list.files(pattern="png"))  
   setwd(wd)  
   system(paste("gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -sOutputFile=",strsplit(pdf,".",fixed=T)[[1]][1],suffix,".pdf ",td,"/*.pdf",sep=""))  
   file.remove(list.files(td,full=T))  
  if(verbose) print("Finished!!")  
 }

MODIS processing with R, GDAL, and the NCO tools

2010-06-03T12:35:00.004-04:00

I use MODIS data for analysis of vegetation dynamics around the world. The native HDF file format provided by NASA is great for archiving the data (it's amazing how much information they include in each file), but unfortunately there aren't many tools for directly (and easily) extracting data stored across files like there are for the NetCDF format. I wanted to be able to use the NCO tools to extract timeseries of a sub-setted region stored in many files. This is quite easy with the ncrcat and similar NCO utilities. The trouble is that the NCO tools require the data to be in netcdf3 or netcdf-4 classic to do this. I tried to use various HDF->netCDF conversion tools but eventually decided it was easiest to convert HDF->geotiff->netcdf. The hierarchical nature of the HDF files makes them difficult to parse and so far the netcdf tools are not able to work with them unless they are 'flattened' to a single group.

So to use the NCO tools with MODIS data, I had to develop a processing routine that will convert them to netcdf-4 classic format. This post will show you how.

I generally use R for my data processing and workflow scripting, though much of this could be done with shell scripting or other languages. The code below is specific to my situation and will almost certainly require some editing if you want to use it. I'm posting it to serve as general guidelines on how one might do this, not as a ready to use and generally applicable set of functions. I appreciate any suggestions on improvements.

First you need to identify which tiles you need by looking at the map:

Then run something like this to download the entire MOD13Q1 (vegetation indices) data from MODIS. Get ready for lots of data (the following list of tiles results in 269GB!):

 ### Some specifications  
 ftpsite="ftp://e4ftl01u.ecs.nasa.gov/MOLT/MOD13Q1.005/" #for veg index  
 hdffolder="1_HDF" #folder to hold HDF files  
 ofolder="2_netcdf" #this is where the output NetCDF files will go  
 bandsubset=paste(1) #keep only these bands from HDF file  
 regionbounds="49 113 30 146" # Set region boundaries for clipping  
 ## Get tiles to download - get these from the map  
 tiles=c("h25v04","h26v04","h26v05","h27v04","h28v04","h27v05","h28v05","h29v05") #list of tiles to download, they will be mosaiced and clipped later - format is important, see wget command below  
 ### download HDFs  
 ### Use wget to recursively download all files for selected tiles - will not overwrite unless source is newer  
 searchstring=paste("\"*",tiles,"*\"",sep="",collapse=",")  
  system(paste("wget -S --recursive --no-parent --no-directories -N -P",hdffolder," --accept",searchstring," --reject \"*xml\"",ftpsite))

This will put all the files in the hdffolder (which I set to "1_HDF"). You now have each tile for time period as a separate file. I find it convenient to mosaic the tiles and perhaps warp them to another projection and subset it to a smaller region than all of the tiles together. To do this I use the following functions (note that some details are specific to extracting NDVI, if you want something else it will need to be edited). The first converts a MODIS date code to the R date format:

modday2date<-function(i) as.Date(paste(c(substr(i,2,5),substr(i,6,9)),collapse="-"),"%Y-%j")
The second converts a geotiff to a netcdf with a single-date record dimension:

 geotiff2netcdf=function(geotiff,output,date=NULL,startday=as.Date("2000-01-01"),varname="NDVI",scale=0.0001,missing=-3000,range=c(-1,1)){ 
  ## Create netcdf file from modistool geotiff output 
  ##create file and open for editing 
  ## convert Date 
  utdates=as.numeric(date-startday) 
  ## Read in geotiff to set netCDF attributes 
  print(paste("Importing ",geotiff)) 
  x=readGDAL(geotiff) 
  ## Set dimentions 
  print("Defining NetCDF") 
  create.nc(output,clobber=T) 
  nc=open.nc(output,write=T) # Opens connection with netcdf file and write permission 
  dim.def.nc(nc, "latitude", dimlength=x@grid@cells.dim[2], unlim=FALSE) 
  dim.def.nc(nc, "longitude", dimlength=x@grid@cells.dim[1], unlim=FALSE) 
  dim.def.nc(nc, "time", dimlength=length(utdates), unlim=TRUE) 
  var.def.nc(nc,"time","NC_SHORT","time") 
  var.put.nc(nc,"time",utdates, start=NA, count=NA, na.mode=0) 
  att.put.nc(nc,"time", "units","NC_CHAR",paste("days since ",startday," 0",sep="")) 
  att.put.nc(nc,"time", "standard_name","NC_CHAR","time") 
  var.def.nc(nc,"longitude","NC_DOUBLE","longitude") 
  att.put.nc(nc,"longitude", "units","NC_CHAR","degrees_east") 
  att.put.nc(nc,"longitude", "standard_name","NC_CHAR","longitude") 
  var.put.nc(nc,"longitude",seq(x@coords[1],x@coords[2],x@grid@cellsize[1]), start=NA, count=NA, na.mode=0) 
  var.def.nc(nc,"latitude","NC_DOUBLE","latitude") 
  att.put.nc(nc,"latitude", "units","NC_CHAR","degrees_north") 
  att.put.nc(nc,"latitude", "standard_name","NC_CHAR","latitude") 
  var.put.nc(nc,"latitude",seq(x@coords[3],x@coords[4],x@grid@cellsize[2]), start=NA, count=NA, na.mode=0) 
  var.def.nc(nc,varname,"NC_SHORT",c("longitude","latitude","time")) 
  att.put.nc(nc,varname, "missing_value","NC_DOUBLE",missing) 
  att.put.nc(nc,varname, "units","NC_CHAR",varname) 
  att.put.nc(nc,varname, "standard_name","NC_CHAR",varname) 
  ## Process the data 
  print("Processing data and writing to disk") 
  notnull=x@data!=missing 
  x@data[notnull][x@data[notnull]range[2]/scale]=missing #get rid of occasional weird value 
  ## write the data to netCDF file 
  var.put.nc(nc,varname,as.matrix(x)[,ncol(as.matrix(x)):1], start=c(1,1,1),count=c(x@grid@cells.dim[1],x@grid@cells.dim[2],1)) 
  close.nc(nc) 
  print("Finished!") 
 }

And the third puts it all together by first using the mrtmosaic tool to mosaic the tiles and subset to the band of interest, then using resample tool (available at the link above) to convert the HDF to a geotiff and then calling the function above to convert that to a netCDF:

 ### Function to procss the files using the commands above.  
 modis2netcdf<-function(i,hdffolder,outputfolder,files,bandsubset=1){  
  tdir="/media/fynbos/public/tmp"  
  dir.create(tdir)  
  tparams=paste(tdir,"/",i,"_params.txt",sep="")  
  tmos=paste(tdir,"/",i,"_mosaic.hdf",sep="")  
  tclip=paste(tdir,"/",i,".tif",sep="")  
  ofile=paste(outputfolder,"/",i,".nc",sep="")  
  f=file.exists(ofile) # file exists?  
  if(f) print(paste(ofile," already exists, moving to next file...."))  # if exists, stop  
  if(!f) {  
  ## Mosaic and subset  
   fs=paste(hdffolder,"/",grep(i,files,value=T),sep="")  
  write.table(fs,file=tparams,row.names=F,col.names=F,quote=F)  
  system(paste("mrtmosaic -i ",tparams," -s ",bandsubset," -o ",tmos,sep=""))  
  ## Clip and reproject  
  if(!file.exists("ModisToolParams.prm")) {print("Parameter file missing, stopping process") ; stop}  
  system(paste("resample -p ModisToolParams.prm -i ",tmos," -o ",tclip," -t GEO -l \'",regionbounds,"\' -a INPUT_LAT_LONG -r CC",sep=""))  
  geotiff=paste(tdir,"/",i,".250m_16_days_NDVI.tif",sep="")  
  geotiff2netcdf(geotiff,output=ofile,date=modday2date(i))  
  file.remove(tmos,tclip,tparams); gc()  
  print(paste(which(dates%in%i), "out of ", length(dates)))  
 }}

You can then run these functions to process a folder of HDF files with the following commands. If you are on a *nix system with multiple cores, I recommend using the multicore package to parallelize the processing. Otherwise you'll need to change the mclapply() below to lapply():

 allfiles=as.character(list.files(path = paste(getwd(),hdffolder,sep="/"))) # build file list  
 dates=unique(do.call("rbind",strsplit(allfiles,".",fixed=T))[,2]) #get unique dates  
 type=do.call("rbind",strsplit(allfiles,".",fixed=T))[1,1] #get prefix of file type  
 band=1  
 ### Run the files - will not overwrite existing mosaics  
 mclapply(dates,modis2netcdf,hdffolder,outputfolder=ofolder,files=allfiles,bandsubset=band,mc.cores=8,mc.preschedule=F)

Now you can use the NCO commands to extract a timeseries from these files, perhaps for only a subset of the data using ncrcat.

Installation of HDF4, HDF5, netCDF4, hdf-eos, and hdf-eos5 on Fedora 11

2010-06-03T11:52:00.001-04:00

I originally wanted to install these libraries to allow processing of MODIS HDF4-EOS files using the NCO operators. I was able to get everything compiled, but the differences between HDF4 and HDF5 and the difficulties of converting the hdf-eos format to the relatively simple (no groups or hierarchical structure) needed for the netcdf-classic format led me to a different processing workflow. More on that later. I'm posting my experience here in the hope that it will save someone the hours I spent figuring this out. The version of netCDF in the Fedora repositories does not seem to have been compiled with NetCDF4/HDF5 support (or szlib encoding) so I decided to compile it all from source. It took a surprising amount of time given that the nco operators depend on netcdf-4, which depends on hdf5 and to convert from hdf4-eos to a version of hdf5-eos that is netcdf compliant requires hdf-eos andhdf-eos5. Here are the steps to get this compiled:

compile szip from source:
./configure --prefix=/usr/local/szip-2.1 --enable-encoding
compile hdf4 from source:
export LD_LIBRARY_PATH=/usr/local/szip-2.1/lib./configure --disable-netcdf --prefix=/usr/local/hdf-4.2.5 --with-szlib=/usr/local/szip-2.1
according to this README, hdf-4 should be compiled without netcdf because it interferes with hdf5 netcdf access.
compile hdf5 from source:
./configure --prefix=/usr/local/hdf5-1.8.4 --with-szlib=/usr/local/szip-2.1
compile netcdf-4 from source with hdf5 support
./configure --prefix=/usr/local/netcdf-4.1.1 --enable-netcdf4 --with-hdf5=/usr/local/hdf5-1.8.4 --with-szlib=/usr/local/szip-2.1 --enable-hdf4=/usr/local/hdf-4.2.5 --with-udunits --enable-hdf4-file-tests --enable-cxx-4 --enable-extra-example-tests CFLAGS=-fPIC
hdf-eos2
./configure CC='/usr/local/hdf-4.2.5/bin/h4cc -Df2cFortran' CFLAGS=-DH5_USE_16_API --prefix=/usr/local/hdfeos2 --with-hdf4=/usr/local/hdf-4.2.5 --enable-install-include --with-szlib=/usr/local/szip-2.1
hdf-eos5
export LD_LIBRARY_PATH=/usr/local/szip-2.1/lib
./configure CC='/usr/local/hdf5-1.8.4/bin/h5cc -Df2cFortran' CFLAGS=-DH5_USE_16_API --prefix=/usr/local/hdfeos5 --with-hdf5=/usr/local/hdf5-1.8.4 --enable-install-include --with-szlib=/usr/local/szip-2.1
compile nco with netcdf4 support (apparently not the default provided in the repository). First check it out of the CVS as explained here
cvs -d:pserver:anonymous@nco.cvs.sf.net:/cvsroot/nco login
cvs -z3 -d:pserver:anonymous@nco.cvs.sf.net:/cvsroot/nco co -r nco-4_0_1 -d nco-4.0.1
Then compile it. The last flag comes from advice found here
./configure --prefix=/usr/local/nco --enable-netcdf-4 HDF5_LIB_DIR=/usr/local/hdf5-1.8.4/lib NETCDF_ROOT=/usr/local/netcdf-4.1.1 NETCDF4_ROOT=/usr/local/netcdf-4.1.1 --disable-dap-netcdf CPPFLAGS='-DHAVE_NETCDF4_H'
Following the instructions at: http://hdfeos.net/software/ncl.php#ref_sec:convert-classic, I wrote a short R program that edits the text file to remove the group information and flatten the file.
compile aug_eos to augment a hdf5 file so it can be accessed with netcdf-4. However this did not allow the NCO tools (which require netcdf-3 or netcdf-4 classic formating).
./configure CPPFLAGS=-DH5_USE_16_API --prefix=/usr/local/aug_eos --with-hdfeos5=/usr/local/hdfeos5 --with-hdf5=/usr/local/hdf5-1.8.4 --with-szlib=/usr/local/szip-2.1