MODIS

Introduction

• About MODIS

There are two MODIS Satellites, Aqua and Terra.

MODIS Aqua

Sea surface color, temperature, and derivative images

• http://oceancolor.gsfc.nasa.gov/PRODUCTS/

The following examples assume you want to import multiple images at once, and so use shell scripting loops.

SST (Level 3)

Get the data

• 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: MODIS/Aqua
* 2002: Year at start
* 353: Julian day at start
* 2002: Year at end
* 360: Julian day at end
* L3m: Level 3 data, mapped (Plate carrée)
* 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)

Import

GRASS can be quite strict about moving maps around between differing map projections. Because GDAL does not automatically georeference HDF4 images, we need to apply the georeferencing information manually. This is a slight chore, but means fewer errors due to projection mis-matches later on. Everything we need to know is given in the data file's meta-data, as seen with gdalinfo.

Method 1

Use GDAL tools to prepare the dataset and create a GeoTiff which can be directly imported into GRASS.

• Create a Lat/Lon GRASS location and mapset

Plate carrée, or WGS84? -- Does it matter with data at the km scale?

-- Encode with +ellps=sphere (?) or epsg:32662?, see http://lists.osgeo.org/pipermail/gdal-dev/2007-September/014134.html or http://thread.gmane.org/gmane.comp.gis.gdal.devel/12666

FIXME: gdalinfo: " Map Projection=Equidistant Cylindrical"

• Use gdal_translate to extract the desired data array from the HDF4 file and convert into a GeoTiff. While creating the GeoTiff specify the map projection, bounding box, and no-data information. Then import this GeoTiff using the r.in.gdal module. In this example we set the projection to be lat/lon WGS84 (EPSG code 4326), even though this may not be entirely correct.

for file in A*SST_4 ; do
  echo "map: $file"

  gdal_translate -a_srs "+init=epsg:4326" -a_nodata 65535 \
    -a_ullr -180 90 180 -90 -co "COMPRESS=PACKBITS" \
    HDF4_SDS:UNKNOWN:"$file":0  ${file}_prep.tif

  r.in.gdal in=${file}_prep.tif out=$file
done

Method 2

Import raw HDF4 data into a XY location, clean it up, then copy into a lat/lon location. This is more work than method 1 but demonstrates using GRASS tools to do the manipulation instead of GDAL tools (the result should be the same).

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, using the appropriate SUBDATASET NAME: '

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

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)

Hack: 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=*SST_4` ; do
  r.null $map setnull=65535
done

Method 3

Like method 2 above, but using i.rectify instead of moving the mapset manually. See the more detailed explanation in the Chlorophyll section below.

Processing

• 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 scaling 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 the 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: MODIS/Aqua
* 2002: Year at start
* 335: Julian day at start
* 2002: Year at end
* 365: Julian day at end
* L3m: Level 3 data, mapped (Plate carrée)
* 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

Method 1

• Create a Lat/Lon GRASS location and mapset

Plate carrée, or WGS84? -- Does it matter with data at the km scale?

• Use gdal_translate to convert the HDF4 data into a GeoTiff, and specify map projection, bounding box, and no-data information. Then import this GeoTiff using the r.in.gdal module. In this example we set the projection to be lat/lon WGS84 (EPSG code 4326), even though this may not be entirely correct.

for file in A*_MO_CHLO_4 ; do
  echo "map: $file"

  gdal_translate -a_srs "+init=epsg:4326" -a_nodata 65535 \
    -a_ullr -180 90 180 -90 -co "COMPRESS=PACKBITS" \
    $file ${file}_prep.tif

  r.in.gdal in=${file}_prep.tif out=$file
done

Method 2

• 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

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)

Hack:: 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=*MO_CHLO_4` ; do
  r.null $map setnull=65535
done

Method 3

Like method 2 above, but instead of moving the mapset directory in the file system use i.rectify to move the maps into the target Lat/Lon location. After importing the images, setting their bounds, and removing NULLs, add all images to an imagery group with the i.group module. Run i.target and select a lat/lon location. The idea here is to run i.rectify with a first-order transform, where the transform uses a scaling factor of 1.0 and rotation of 0.0, i.e. no change at all. The trick is to set those. Usually you would use i.points or the GUI georeferencing tool to set those, and here you can add a few points by hand if you like. But the POINTS file in the group's data directory should use identical coordinates for both source and destination. Hack: edit the POINTS file by hand to make it so, using the four corners and 0,0 at the center of the image. Finally run i.rectify to push the image into the target location.

Processing

• 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

• Check:

r.info A20023352002365.L3m_MO_CHLO_4.chlor_A
r.univar A20023352002365.L3m_MO_CHLO_4.chlor_A

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 mg/m^3.

• 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

SeaWiFS

From OceanColor gsfc.nasa.gov

• the SeaWiFS Standard Mapped Products are already geocoded in WGS84 so should be easy to load into a lat/lon location with r.in.gdal. (at least this is true for PAR maps)
• Processing is very similar to the MODIS examples above.

Pathfinder SST

About

AVHRR sea surface temperature data v5.0 at 4km resolution is available for 1985-present.

• See the 4km AVHRR Pathfinder homepage at NOAA's NODC. (especially the user's guide)

• The data can be downloaded from their FTP site:

ftp://data.nodc.noaa.gov/pub/data.nodc/pathfinder/

• The filesize for a 4km resolution data file is approx 16mb.

• You will probably want to download both the raw sst data and the quality mask layer ("*-sst" and "*-qual" files). During processing we will use the quality layer to mask the raw data for valid pixels.

File names look like: 2000301-2000305.s0451pfv50-sst-16b.hdf

key:

* 2000: Year at start
* 301: Julian day at start
* 2000: Year at end
* 305: Julian day at end
* s: 16-bit data
* 04: 4km resolution (8192x4096 image, approx 2.637 minute resolution)
* 5: 5 day ensemble
* 1: nighttime (3 is daytime)
* pfv50: Pathfinder version 5.0
* sst: Product: Raw SST (including cloud hits)
* 16b: 16 bits (only present in pre-2003 data)

Check

• Use GDAL's gdalinfo tool to view the HDF file's meta-data, as above for MODIS HDF4 data.

Import

• This data is already in Lat/Lon WGS84, so things are made a bit easier.

• Open a lat/lon WGS84 location (EPSG code 4326)

• Use gdal_translate to convert the HDF4 file into a GeoTiff containing map projection, bounding box, and no-data information. Then import this GeoTiff using the r.in.gdal module.

# Convert and import SST maps
for file in *pfv50-sst*.hdf ; do
  echo "map: $file"

  gdal_translate -a_srs "+init=epsg:4326" -a_nodata 0 \
    -a_ullr -180 90 180 -90 -co "COMPRESS=PACKBITS" \
    HDF4_SDS:UNKNOWN:"$file":0  ${file}_prep.tif

  r.in.gdal in=${file}_prep.tif out=`basename $file .hdf`
done


# Convert and import Overall Quality layers
for file in *pfv50-qual.hdf ; do
  echo "map: $file"

 gdal_translate -a_srs "+init=epsg:4326" \
   -a_ullr -180 90 180 -90 -co "COMPRESS=PACKBITS" \
   HDF4_SDS:UNKNOWN:"$file":0  ${file}_prep.tif

  map=`basename $file .hdf`
  r.in.gdal in=${file}_prep.tif out="$map"
  r.support "$map" title="SST Overall Quality" history=""
  r.support "$map" history="key: 0 is worst, 7 is best"
  r.colors -r "$map"  # use standard colors
done

Processing

• Convert values to temperature in degrees C. Note slope and intercept are taken from the image's metadata using gdalinfo:

 scale_factor=0.075
 add_off=-3

Given those, we can now process the data. Note the special care is taken to quote the map names in r.mapcalc as they contain minus signs.

for map in `g.mlist rast pat=*pfv50-sst*` ; do
  echo "$map"
  g.region rast=$map
  #r.info -r $map
  Scale=0.075
  Intercept=-3
  r.mapcalc "\"${map}.degC\" = ($Scale * \"$map\") + $Intercept"
  r.support "${map}.degC" units="deg-C"
  #r.colors "${map}.degC" color=bcyr  # fairly nice
done

• Check:

r.info 2000301-2000305.s0451pfv50-sst-16b.degC
r.univar 2000301-2000305.s0451pfv50-sst-16b.degC

• With r.univar (or d.histogram) we see some spurious data points (125deg C water!). We can check with r.univar's "-e percentile=99" that the majority of the data really does not go above 30deg C. If you like you can use r.mapcalc to remove anything above 40C:

r.mapcalc "clean_map = if( orig_map > 40, null(), orig_map )"

Set colors

For a very nice looking image we can use Goddard's OceanColor palettes from

http://oceancolor.gsfc.nasa.gov/PRODUCTS/colorbars.html

For more information on how to use those with GRASS see the above section regarding MODIS data. To apply the converted rules: (download them here)

r.colors "${map}.degC" rules=palette_sst.gcolors

Removing holes

Cloud cover or technical problems often causes there to be some gaps in the data. GRASS has a r.fillnulls module which uses spline interpolation to fill in these small gaps. These holes are left in the data for a reason, so check if it is appropriate for your analysis before removing them. Also note that the areas around the holes may be of low quality and using them as the starting point for filling the holes just spreads the error wider.

The first step is to create a mask of the sea, i.e. the area which it is valid to interpolate over. To do this we overlay a time series of satellite images and see which areas are permanently no data.

• Overlay all maps; the data will be mixed up and bogus but the spatial coverage will be informative

r.patch in=`g.mlist -r rast pat='A200.*4$' sep=,` out=all_sat_images_stacked

• Create a 0/1 mask layer showing areas of sea, and label categories

r.mapcalc "sea_mask = if( isnull(all_sat_images_stacked), 0, 1)"
r.category sea_mask rules=- << EOF
  0:Land
  1:Sea
EOF

g.remove all_sat_images_stacked

The next step is to perform the filling operation. To speed up processing, we first temporarily fill known land areas with a value as there is no point interpolating over those gaps. To avoid interpolation overshoots, smear coastal values slightly inland before interpolation. After the holes are filled we reapply the mask to again remove those areas. As the r.fillnulls step can take some time, it is advisable to use g.region to zoom into an area of interest first.

for map in `g.mlist rast pat=*SST_4.degC` ; do
   r.grow in=$map out=${map}.grown radius=5
   r.mapcalc "tmp_landfill = if( isnull($map) && sea_mask == 1, \
     null(), if( isnull(${map}.grown), 0, ${map}.grown) )"
   r.fillnulls in=tmp_landfill out=${map}.filled_tmp
   g.copy sea_mask,MASK
   r.mapcalc "${map}.filled = ${map}.filled_tmp"
   g.remove tmp_landfill,${map}.grown,${map}.filled_tmp,MASK
done