MODIS

From GRASS-Wiki
Jump to navigation Jump to search

There are two MODIS Satellites, Aqua and Terra.

MODIS Aqua

Sea surface color, temperature, and derivative images

SST (Level 3)

Get data:

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: 4.6km pixel size (8640x4320 image, 2.5 minute resolution)


Decompress: 

bzip2 -d A20023532002360.L3m_8D_SST_4.bz2

Check:

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

Example of SST metadata from gdalinfo: 

 Parameter=Sea Surface Temperature
 Measure=Mean
 Units=deg-C
 Scaling Equation=(Slope*l3m_data) + Intercept = Parameter value
 Slope=0.000717185
 Intercept=-2
 Scaled Data Minimum=-2
 Scaled Data Maximum=45
 Data Minimum=-1.999999
 Data Maximum=37.06

Subdatasets: 

 SUBDATASET_1_NAME=HDF4_SDS:UNKNOWN:"A20023532002360.L3m_8D_SST_4":0
 SUBDATASET_1_DESC=[4320x8640] l3m_data (16-bit unsigned integer)
 SUBDATASET_2_NAME=HDF4_SDS:UNKNOWN:"A20023532002360.L3m_8D_SST_4":1
 SUBDATASET_2_DESC=[4320x8640] l3m_qual (8-bit unsigned integer)


Extract:

Because the HDF4 data file contains multiple data arrays we need to extract the one we are interested in. In this example we extract the data layer but not the quality layer, and output a GeoTIFF file using GDAL's gdal_translate tool: 

for file in *SST_4 ; do
  gdal_translate HDF4_SDS:UNKNOWN:"$file":0 ${file}_unproj.tif
done


Import:

  • Create a simple XY GRASS location and mapset
  • Import the imagery and set image bounds and resolution:
for file in *unproj.tif ; do
  BASE=`basename $file _unproj.tif`
  echo "map: $file"
  r.in.gdal in=$file out=$BASE
  r.region $BASE n=90 s=-90 w=-180 e=180
done
  • Check:
g.list rast
r.info A20023532002360.L3m_8D_SST_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 temperature in degrees C. Note slope and intercept are taken from the image's metadata using gdalinfo.
for map in `g.mlist rast pat=*SST_4` ; do
  echo "$map"
  g.region rast=$map
  #r.info -r $map
  Slope=0.000717185
  Intercept=-2
  r.mapcalc "${map}.degC = ($Slope * $map) + $Intercept"
  r.support "${map}.degC" units="deg-C"
  #r.colors "${map}.degC" color=bcyr  # fairly nice
done


  • Check:
r.info A20023532002360.L3m_8D_SST_4.degC
r.univar A20023532002360.L3m_8D_SST_4.degC

Set colors:

Although the standard bcyr color map looks nice, for a very nice image 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 temperatures are -2-45 deg C.

  • 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 temperature values
echo "# Color rules for MODIS SST" > palette_sst.gcolors
Slope=0.000717185
Intercept=-2
cat palette_sst.txt | grep -v '^#' | \
  awk -v Slope=$Slope -v Intercept=$Intercept \
    '{ printf("%f %d:%d:%d\n", \
        (Slope * (($1 +1)^2 -1) + Intercept), $2, $3, $4)}' \
  >> palette_sst.gcolors
echo "nv black" >> palette_sst.gcolors
# better: edit last rule to be 45.000719 206:206:206


The processed color rules file is here: [1]


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

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 are logarithmic which poses some challenges to rendering nice colors. (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: [2]


  • 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
Retrieved from "https://grasswiki.osgeo.org/w/index.php?title=MODIS&oldid=6139"