MODIS

2009-05-17T21:40:29Z

⚠️Arctos: /* MODIS Aqua */ Specified input file assumption

== Introduction ==
* [http://modis.gsfc.nasa.gov/ About MODIS]

There are two satellites, Aqua and Terra which carry the MODIS sensor as payload. The Moderate Resolution Imaging Spectroradiometer (MODIS) is a 36-channel from visible to thermal-infrared sensor that was launched as part of the Terra satellite payload in December 1999 and Aqua satellite (May 2002). The Terra satellite passes twice a day (at about 10:30am, and 22:30pm local time), also the Aqua satellite passes twice a day (at about 01:30am, and 13:30pm local time).

MODIS Sinusoidal Grid SHAPE files can be downloaded [http://gis.fem-environment.eu/modis-sinusoidal-gis-files/ here]

== MODIS Terra ==

MODIS on Terra provides many data products including:
* Daily Surface Reflectance
* Snow coverage
* NDVI / EVI Vegetation index
* Land Surface Temperature (LST)
* Atmosphere Daily Global Products
* many more

See [http://modis-land.gsfc.nasa.gov/ MODIS Land Team Home Page]

=== Advanced MODIS LST time series reconstruction ===

* See http://gis.fem-environment.eu/modis-lst

== 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. They also assume input files are SMI (Standard Mapped Image) format, as produced by the Ocean Color group (http://oceancolor.gsfc.nasa.gov).

=== 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: <tt>A20023532002360.L3m_8D_SST_4</tt>
:''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 [http://www.gdal.org GDAL's] <tt>gdalinfo</tt> 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 <tt>gdalinfo</tt>.

Alternatively the USGS/NASA's [https://lpdaac.usgs.gov/lpdaac/tools/modis_reprojection_tool MODIS Reprojection Tool] (MRT) for Level 3 data or [https://lpdaac.usgs.gov/lpdaac/tools/modis_reprojection_tool_swath MODIS Swath Reprojection Tool] for Level 1B or 2 data can be used to create a GeoTiff.

===== 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 "<tt>+nadgrids=@null</tt>" trick?? http://proj.maptools.org/faq.html#sphere_as_wgs84

* Use <tt>gdal_translate</tt> 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 <tt>r.in.gdal</tt> 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 <tt>gdal_translate</tt> 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 <tt>proj:</tt> 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 <tt>i.rectify</tt> 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 <tt>bcyr</tt> 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: [http://trac.osgeo.org/grass/browser/grass-addons/raster/r.colors.tools/palettes]

* 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: <tt>A20023352002365.L3m_MO_CHLO_4</tt>
:''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 [http://www.gdal.org GDAL's] <tt>gdalinfo</tt> 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 <tt>gdal_translate</tt> to convert the HDF4 data into a GeoTiff, and specify map projection, bounding box, and no-data information. Then import this GeoTiff using the <tt>r.in.gdal</tt> 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 <tt>proj:</tt> 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 <tt>i.group</tt> module. Run <tt>i.target</tt> 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 <tt>i.points</tt> 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: [http://trac.osgeo.org/grass/browser/grass-addons/raster/r.colors.tools/palettes]

* 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

=== Removing holes ===

See [[AVHRR]].

=== Rolling composite ===

Create a rolling 28 day composite:

Say you have weekly satellite data but there is a lot of NULL cells where there was cloud cover that day. You can fill these holes by looking for valid data in the previous week's image, and if again you find clouds, look in week prior to that, etc., up to 4 weeks have past, after which point if any holes remain give up and leave it unfilled. You can do this quite simply with r.patch:

r.patch output=rolling28day.date input=newest,last_week,2weeks_ago,3weeks_ago

NULL data cells are filled using the maps in the given '''input=''' order.

== See also ==

* [[SeaWiFS]]
* [[AVHRR]]

Map Reprojection

2007-03-15T21:04:45Z

⚠️Arctos:

'''Q:''' How to change map/image projections, datums, etc in GRASS GIS?

'''A:'''

# create an import location in the projection and datum of the map/image you want to import
# import the map/image into the location
# create a destination location in the projection and datum you want to reproject the map/image into
# working in the destination location, use r.proj or v.proj (depending on whether the map/image is raster or vector) to reproject the map/image from the import location to the destination location.

If your raster data set is too big (i.e, 2 gb or so), r.proj will crash. Instead, export as a geotif, use gdalwarp to change the project and then import the new geotif file. If you have installed gdal, you almost certainly have gdalwarp as well. Directions for gdalwarp are at http://www.remotesensing.org/gdal/gdalwarp.html. The options for gdalwarp are a bit confusing. Here is an example of projecting a geotif based on srtm elevation data to utm 37N. -t_srs is the output file projection. AfricaHornElev.tif is the input file, AfricaHornElev37n.tif is the output file.

gdalwarp -t_srs '+proj=utm +zone=37 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' -rcs AfricaHornElev.tif AfricaHornElev37n.tif

Another (similar) strategy frequently needed when maps need to be imported using, e.g. v.in.ogr or r.in.ogr, is to tell these programs to create a new location (argument "location"). If the imported data contains proper PROJ_INFO files, they can be reprojected into the current mapset, using the the projection in the current location/mapset.
[[Category:FAQ]]

GRASS-Wiki - User contributions [en]

MODIS

Map Reprojection