From GRASS-Wiki
Revision as of 03:34, 1 April 2008 by HamishBowman (talk | contribs) (Chlorophyll-a (Level 3))

Jump to: navigation, search

There are two MODIS Satellites, Aqua and Terra.


Sea surface color, temperature, and derivative images

SST (Level 3)

File names look like: A20023532002360.L3m_8D_SST_4 key:

* A: Aqua
* 2002: Year at start
* 353: Julian day at start
* 2002: Year at end
* 360: Julian day at end
* L3m: Level 3 data, mapped
* 8D: 8 day ensemble
* SST: Sea Surface Temperature product
* 4: 4km pixel size (8640x4320 image, 2.5 minute resolution)

Chlorophyll-a (Level 3)

Get data:

File names look like: A20023352002365.L3m_MO_CHLO_4

* A: Aqua
* 2002: Year at start
* 335: Julian day at start
* 2002: Year at end
* 365: Julian day at end
* L3m: Level 3 data, mapped
* MO: One month ensemble
* CHLO: Chlorophyll a concentration product
* 4: 4.6km pixel size (8640x4320 image, 2.5 minute resolution)


bzip2 -d A20023352002365.L3m_MO_CHLO_4.bz2


  • Use GDAL's gdalinfo tool to view the HDF file's meta-data.
gdalinfo A20023352002365.L3m_MO_CHLO_4

Example of CHLO metadata from gdalinfo:

 Parameter=Chlorophyll a concentration
 Units=mg m^-3
 Scaling Equation=Base**((Slope*l3m_data) + Intercept) = Parameter value
 Scaled Data Minimum=0.01
 Scaled Data Maximum=64.5654
 Data Minimum=0.002637
 Data Maximum=99.99774


  • Create a simple XY GRASS location and mapset
  • Import the imagery and set image bounds and resolution:
for file in A*_MO_CHLO_4 ; do
  echo "map: $file" in=$file out=$file
  r.region $file n=90 s=-90 w=-180 e=180
  • Check:
g.list rast A20023352002365.L3m_MO_CHLO_4
  • Convert the simple XY location into a Lat/Lon location
(work in progress: we badly need a cleaner method)
Integerized Sinusoidal (ISIN) projection, Plate carrée, or WGS84? -- Does it matter with data at the km scale?
Modify the location's projection settings with g.setproj?
What I did: after import move the mapset dir into a world_ll location, then edit the mapset's WIND file and change the proj: line from 0 to 3. (0 is unprojected, 3 is LL) You also have to edit each file header in the $MAPSET/cellhd/ directory.
  • Remove NULLs
for map in `g.mlist rast pat=A200*` ; do
  r.null $map setnull=65535
  • Convert to chlorophyll a concentration. Note slope and intercept are taken from the image's metadata using gdalinfo.
for map in `g.mlist rast pat=*MO_CHLO_4` ; do
  echo "$map" -r $map
  r.mapcalc "${map}.chlor_A = 10^(($Slope * $map) + $Intercept)" "${map}.chlor_A" units="mg m^-3"
  #r.colors -e "${map}.chlor_A" color=bcyr  # :-(

Set colors:

The Chlorophyll maps a logarithmic which poses some challenges to rendering it nicely. (the unconverted import image displays nicely with a linear color map) We can use Goddard's OceanColor palettes from

Those are given in 0-255 range, the imported HDF4 data is 0-65535, and the converted chlorophyll a concentration is 0-65.

  • Use the UNIX awk tool to convert each color rule with the same exponential function we used to convert the data, and save it to a rules file:
# scale 0-255 to 0-65535 and then convert to chlor-a values
echo "# Color rules for MODIS Chloropyll-a" > palette_chl_etc.gcolors
cat palette_chl_etc.txt | grep -v '^#' | \
  awk -v Slope=$Slope -v Intercept=$Intercept \
    '{ printf("%f %d:%d:%d\n", \
        10^((Slope * (($1 +1)^2 -1)) + Intercept), $2, $3, $4)}' \
  >> palette_chl_etc.gcolors
echo "nv black" >> palette_chl_etc.gcolors
# better: edit last rule to be 64.574061 100:0:0

The processed color rules file is here: [1]

  • Apply color rules to imported maps:
for map in `g.mlist rast pat=*MO_CHLO_4` ; do
  r.colors "${map}.chlor_A" rules=palette_chl_etc.gcolors