Plotting MODIS data in Python / basemap - a MODIS workflow solution?

[John [H2O]]

Hello,

Recently I read a clear and helpful blog entry by Christian Perone (author of pyevolve):
http://pyevolve.sourceforge.net/wordpress/?p=86

I asked him if he had a solution for plotting standard MODIS hdf products available here:
http://rapidfire.sci.gsfc.nasa.gov/realtime/2009300/

In more general terms, what are people using for a 'MODIS workflow'? I can't imagine I'm the first to want to plot MODIS images using basemap. Does anyone have a 'roadmap' (as Christian phrased it)? I am interested in the steps from A-Z. That is, what tool to read the hdf files (modis is hdf4), what tool to reproject the data, and finally, the basemap plotting.

Looking forward to responses,
john

[Vincent Schut-2]

John [H2O] wrote:

> Hello,
>
> Recently I read a clear and helpful blog entry by Christian Perone (author
> of pyevolve):
> http://pyevolve.sourceforge.net/wordpress/?p=86
>
> I asked him if he had a solution for plotting standard MODIS hdf products
> available here:
> http://rapidfire.sci.gsfc.nasa.gov/realtime/2009300/
>
> In more general terms, what are people using for a 'MODIS workflow'? I can't
> imagine I'm the first to want to plot MODIS images using basemap. Does
> anyone have a 'roadmap' (as Christian phrased it)? I am interested in the
> steps from A-Z. That is, what tool to read the hdf files (modis is hdf4),
> what tool to reproject the data, and finally, the basemap plotting.

FWIW, I'm using gdal for both reading and reprojecting (working with
modis daily spectral images, mainly, e.g. MOD09GA/Q, in sinusoidal
projection). I implement everything in python (gdal python bindings,
numpy, scipy, etc), including processing, so I don't use the default
gdal utilities much for this. Also I don't plot, I just output geotiffs
which I view using openev or GIS software.

Vincent.

>
> Looking forward to responses,
> john
>
>
>


------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

[Jose Gomez-Dans]

Hi,

I asked him if he had a solution for plotting standard MODIS hdf products
available here:
http://rapidfire.sci.gsfc.nasa.gov/realtime/2009300/

In more general terms, what are people using for a 'MODIS workflow'? I can't
imagine I'm the first to want to plot MODIS images using basemap. Does
anyone have a 'roadmap' (as Christian phrased it)? I am interested in the
steps from A-Z. That is, what tool to read the hdf files (modis is hdf4),
what tool to reproject the data, and finally, the basemap plotting.

Really, it's what you want to do with your MODIS data. My "workflow" is usually as follows:
1.- Access MODIS data (and ancillary stuff, such as QA flags etc) using Python's GDAL bindings.
2.- Manipulate the MODIS data from (1) using numpy, scipy. If there's significant looping involved, use weave to speed things up.
3.- Plot using matplotlib. Usually, as imshow ("vanilla matplotlib"), sometimes using basemap. The difference is whether I'm just quickly plotting something together, or whether I want to actually have a map where I want to plot other stuff on top of the MODIS data.

In terms of rapidfire, you can even access the WMS data using GDAL, so no need to download the data. An example for monitoring El Niño related fires in Borneo is attached. You need to get the fires in the last 7days file (we download it as a crontab jobby), and the attached script just plots the data on a basemap for quick visualisation.

I am in the process of putting most of my notes on <http://sites.google.com/site/spatialpython/> and there's also the Unofficial Python GIS SIG <http://groups.google.com/group/python-gis-sig>, another useful resource.

Your particular example is plotting RGB composites derived from Level 1 data, it seems? As I said, if you tell us what you want to do with it, we may be able to provide more information.

Cheers,

Jose

[do_borneo.py]

---

# -*- coding: utf-8 -*-
#Copyright: J Gomez-Dans <j.gomez-dans@...>
#
#2009-10-27
#
#This code is licensed under the terms of the GNU General Public License (GPL),
#version 2 or above; See /usr/share/common-licenses/GPL , or see
#http://www.fsf.org/licensing/licenses/gpl.html

from osgeo import ogr
import pylab
from mpl_toolkits.basemap import Basemap

g = ogr.Open("SouthEast_Asia_7d.shp")
L = g.GetLayer(0)
x=[] ; y=[]
while True:
        feat = L.GetNextFeature()
        if not feat: break
        geom = feat.GetGeometryRef()
        x.append( geom.GetX()) ; y.append( geom.GetY())
m = Basemap( llcrnrlon=108, llcrnrlat=-5, urcrnrlon=120,urcrnrlat=8, projection='merc',\
                        lat_ts=0.0,resolution='h')
m.bluemarble()
m.drawcoastlines (linewidth=0.25,color='w')
m.drawcountries(linewidth=0.25,color='w')
m.drawmeridians(pylab.arange(0,360,5))
m.drawparallels(pylab.arange(-90,90,5))
m.drawmapboundary()
(X,Y) = m( x, y )

m.plot ( X, Y , 'rs', ms=7.0, mew=0)
pylab.title("Fires in Borneo in the last 7 days")
pylab.savefig ("Borneo.png")

---

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

[John [H2O]]

Jose Gómez-Dans-2 wrote:

Really, it's what you want to do with your MODIS data. My "workflow" is
usually as follows:
1.- Access MODIS data (and ancillary stuff, such as QA flags etc) using
Python's GDAL bindings.
2.- Manipulate the MODIS data from (1) using numpy, scipy. If there's
significant looping involved, use weave to speed things up.
3.- Plot using matplotlib. Usually, as imshow ("vanilla matplotlib"),
sometimes using basemap. The difference is whether I'm just quickly plotting
something together, or whether I want to actually have a map where I want to
plot other stuff on top of the MODIS data.


Your particular example is plotting RGB composites derived from Level 1
data, it seems? As I said, if you tell us what you want to do with it, we
may be able to provide more information.

Cheers,
Jose

Thanks Jose and Vincent for the quick replies. I'll have to read up on the gdal bindings it seems.

What I want to do is quite simple (I think). I simply would like to be able to:

1) Download any resolution of the hdf files available via :
http://rapidfire.sci.gsfc.nasa.gov/realtime/2009XXX/  where XXX is day of year.

2) Import the hdf into python, access the data (reproject if necessary). Yes, it's the merged RGB level 1 data I want to ultimately plot.

2b) produce geotiffs for use elsewhere

3) Plot the data in basemap, plot some other data (my own, xy tracks, points, etc) over the data.

4) Save png files


I'll have a look at your notes now.

-john

[Jose Gomez-Dans]

Hi!

This is probably getting a bit off-topic now, maybe better keep this converation going on python-gis-sig list? I'm cc'ing there....

2009/10/27 John [H2O] <washakie@...>

What I want to do is quite simple (I think). I simply would like to be able
to:

1) Download any resolution of the hdf files available via :
http://rapidfire.sci.gsfc.nasa.gov/realtime/2009XXX/  where XXX is day of
year.

Basically you have the MOD021K, MOD02HK and MOD02QK data granules, plus the MOD03 data granule. The former are the 1k, 500m and 250m data (hence the name) and the latter is a collection of geolocation terms (imaging geometries, etc).

You realise that these are swath data don't you? So they have no direct information on the location of each pixel, although you can work it out from the sensor height, attitude and so on. 

2) Import the hdf into python, access the data (reproject if necessary).
Yes, it's the merged RGB level 1 data I want to ultimately plot.

Importing into python is easy with GDAL. Your reprojection is a problem because the data aren't even projected. As I said above, they are swath observations that need to be gridded (Level 2 product).

2b) produce geotiffs for use elsewhere

Again, easy using GDAL (provided you don't care about using swath data, OR that you know how to grid it).

3) Plot the data in basemap, plot some other data (my own, xy tracks,
points, etc) over the data.

Again, you need to know the projection of the MODIS data for this to be of any use.

I don't have an off-hand way of gridding the Level 1 data, but someone may have. Other than the MS2GT <http://nsidc.org/data/modis/ms2gt/> (MODIS Swath to Grid toolbox), there's also the information given in <http://www.mcst.ssai.biz/mcstweb/L1B/product.html>, and I guess you could use that as a way of developing a python MS2GT ;D

Jose

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users