MODIS
There are two MODIS Satellites, Aqua and Terra.
MODIS Aqua
Sea surface color, temperature, and derivative images
SST (Level 3)
- Get data in HDF4 format from http://oceancolor.gsfc.nasa.gov/PRODUCTS/L3_sst.html
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:
- Get data in HDF4 format from http://oceancolor.gsfc.nasa.gov/PRODUCTS/L3_chlo.html
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