MODIS: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
Line 91: Line 91:
  for map in `g.mlist rast pat=*MO_CHLO_4` ; do
  for map in `g.mlist rast pat=*MO_CHLO_4` ; do
   echo "$map"
   echo "$map"
  g.region rast=$map
   #r.info -r $map
   #r.info -r $map
   Slope=5.813776e-05
   Slope=5.813776e-05

Revision as of 10:36, 1 April 2008

There are two MODIS Satellites, Aqua and Terra.

MODIS Aqua

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

key:
* 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)


Decompress:

bzip2 -d A20023352002365.L3m_MO_CHLO_4.bz2


Check:

  • 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
 Measure=Mean
 Units=mg m^-3
 Scaling=logarithmic
 Scaling Equation=Base**((Slope*l3m_data) + Intercept) = Parameter value
 Base=10
 Slope=5.813776e-05
 Intercept=-2
 Scaled Data Minimum=0.01
 Scaled Data Maximum=64.5654
 Data Minimum=0.002637
 Data Maximum=99.99774

Import:

  • 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"
  r.in.gdal in=$file out=$file
  r.region $file n=90 s=-90 w=-180 e=180
done
  • Check:
g.list rast
r.info 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? (run that from the PERMANENT mapset)

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
done
  • 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"
  g.region rast=$map
  #r.info -r $map
  Slope=5.813776e-05
  Intercept=-2
  r.mapcalc "${map}.chlor_A = 10^(($Slope * $map) + $Intercept)"
  r.support "${map}.chlor_A" units="mg m^-3"
  #r.colors -e "${map}.chlor_A" color=bcyr  # :-(
done


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 http://oceancolor.gsfc.nasa.gov/PRODUCTS/colorbars.html

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
Slope=5.813776e-05
Intercept=-2
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
done