[GRASS-user] GRASS & MODIS L3 data

GRASS GRASS user
1

  [Since my posting seem not to be delivered, I'm resending it.
  As a subscribed person, this time.]

  Browsing through the list, I've got an impression that GRASS
  cannot import MODIS Level 2 data directly. IIUC, the main
  problem is that the consequent rows of the data sets overlap
  ``on the Earth'', and it couldn't be fixed by either `i.rectify'
  or `gdalwarp'. Please correct me if I'm wrong here.

  Is the use of a reprojection tool the only solution? I know
  about MODIS Reprojection Tool [1] and MODIS Swath-to-Grid
  Toolbox [2]; are there any other tools available? Have anyone
  made a comparison between these tools?

[1] http://edcdaac.usgs.gov/landdaac/tools/modis/index.asp
[2] http://nsidc.colorado.edu/data/modis/ms2gt/

2

For myself, I just use HEGTools which is freely available from NASA and
convert everything to GeoTiff files. It allows you to select a variety of
projections, create multiband composites, or produce single bands. You
can also sub-sample the data set.

Best, MEH

3

Hi Ivan,

On 30 Mar 2007 22:24:36 +0700, Ivan Shmakov <ivan@ivanshmakov.homeip.net> wrote:

        Browsing through the list, I've got an impression that GRASS
        cannot import MODIS Level 2 data directly. IIUC, the main
        problem is that the consequent rows of the data sets overlap
        ``on the Earth'', and it couldn't be fixed by either `i.rectify'
        or `gdalwarp'. Please correct me if I'm wrong here.

I use a lot of MODIS data (all levels essentially). Rather than using
the HEGTools, I have written a small python script that uses GDAL (and
Markus Neteler's notes on how to reproject MODIS data) to reproject
and select the area I want. The ouput is GeoTIFF, which can then be
easily imported into GRASS.

If you want to have a look at it, let me know, and I'll send it over.

Jose

4

Jose Gomez-Dans <jgomezdans@gmail.com> writes:

>> Browsing through the list, I've got an impression that GRASS cannot
>> import MODIS Level 2 data directly. IIUC, the main problem is that
>> the consequent rows of the data sets overlap ``on the Earth'', and
>> it couldn't be fixed by either `i.rectify' or `gdalwarp'. Please
>> correct me if I'm wrong here.

> I use a lot of MODIS data (all levels essentially). Rather than using
> the HEGTools, I have written a small python script that uses GDAL
> (and Markus Neteler's notes on how to reproject MODIS data)

  Where I could find these notes?

> to reproject and select the area I want. The ouput is GeoTIFF, which
> can then be easily imported into GRASS.

> If you want to have a look at it, let me know, and I'll send it over.

  It would be very interesting to me to see it.

5

"JG" == Jose Gomez-Dans <jgomezdans@gmail.com> writes:

  Looks like I've mistyped the Subject:. Actually, I'm interested
  in *both* Level 2 (scan-wise) and Level 3 (sinusoidal
  projection) MODIS data.

[...]

> I use a lot of MODIS data (all levels essentially). Rather than using
> the HEGTools, I have written a small python script that uses GDAL
> (and Markus Neteler's notes on how to reproject MODIS data)

  I. e.:

http://mpa.itc.it/markus/useful/modis_hdf2erdas_ll_wgs84.sh

  I seem not to understand the hack (why +ellps=sphere in
  -t_srs?), but it, indeed, works. Thanks.

> to reproject and select the area I want.

[...]

  Although `proj -b' and MS2GT's `fornav' work reasonably fine,
  I'm still interested in any other tools that might help in
  importing MODIS L2 data.

6

Ivan Shmakov wrote:

"JG" == Jose Gomez-Dans <jgomezdans@gmail.com> writes:

  Looks like I've mistyped the Subject:. Actually, I'm interested
  in *both* Level 2 (scan-wise) and Level 3 (sinusoidal
  projection) MODIS data.

[...]

> I use a lot of MODIS data (all levels essentially). Rather than using
> the HEGTools, I have written a small python script that uses GDAL
> (and Markus Neteler's notes on how to reproject MODIS data)

  I. e.:

http://mpa.itc.it/markus/useful/modis_hdf2erdas_ll_wgs84.sh

  I seem not to understand the hack (why +ellps=sphere in
  -t_srs?), but it, indeed, works.

This is related I guess:
http://proj.maptools.org/faq.html#sphere_as_wgs84

Maciek

7

Maciej Sieczka wrote on 04/03/2007 06:57 PM:

Ivan Shmakov wrote:
  

"JG" == Jose Gomez-Dans <jgomezdans@gmail.com> writes:
              

  Looks like I've mistyped the Subject:. Actually, I'm interested
  in *both* Level 2 (scan-wise) and Level 3 (sinusoidal
  projection) MODIS data.

[...]

> I use a lot of MODIS data (all levels essentially). Rather than using
> the HEGTools, I have written a small python script that uses GDAL
> (and Markus Neteler's notes on how to reproject MODIS data)

  I. e.:

http://mpa.itc.it/markus/useful/modis_hdf2erdas_ll_wgs84.sh

  I seem not to understand the hack (why +ellps=sphere in
  -t_srs?), but it, indeed, works.
    
This is related I guess:
http://proj.maptools.org/faq.html#sphere_as_wgs84
  

I don't think that it is related.
I think that it is a GDAL issue. You may check and join my bug report
http://trac.osgeo.org/gdal/ticket/1239

Since the hack mostly works, I don't bother any more :slight_smile:
In most cases I am also using the MODIS Reprojection Tool.

Markus

------------------
ITC -> dall'1 marzo 2007 Fondazione Bruno Kessler
ITC -> since 1 March 2007 Fondazione Bruno Kessler
------------------

8

"MN" == Markus Neteler <neteler@itc.it> writes:

[...]

>>> I seem not to understand the hack (why +ellps=sphere in -t_srs?),
>>> but it, indeed, works.

>> This is related I guess:
>> http://proj.maptools.org/faq.html#sphere_as_wgs84

> I don't think that it is related. I think that it is a GDAL
> issue. You may check and join my bug report
> http://trac.osgeo.org/gdal/ticket/1239

[...]

  Well, I've first encountered the ``MODIS L3'' problem when I've
  tried to import several files into GRASS using `r.in.gdal'.
  After a few vector maps were reprojected into the location as
  well, it was noticed, that either rasters or vectors were
  shifted (or both.)

  Now, I think that it's an issue with `v.proj', while `r.in.gdal'
  actually works fine. Please consider the following (tested on
  GRASS 6.0.2 debian 6.)

  Let me introduce two locations first:

world-ll $ g.proj -p
-PROJ_INFO-------------------------------------------------
name : Latitude-Longitude
datum : wgs84
towgs84 : 0.000,0.000,0.000
proj : ll
ellps : wgs84
-PROJ_UNITS------------------------------------------------
unit : degree
units : degrees
meters : 1.0
world-ll $

modis-sinu $ g.proj -p
-PROJ_INFO-------------------------------------------------
name : Sinusoidal (Sanson-Flamsteed)
towgs : 0,0,0
proj : sinu
R : 6371007.181
a : 6371007.181
b : 6371007.181
es : 0.0000000000
f : 0.0000000000
lat_0 : 0.0000000000
lon_0 : 0.0000000000
-PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1.0
modis-sinu $

  ... And a `points' file (longitude-latitude coordinates):

world-ll $ cat points-ll.dat
80.00|50.00|a point on L3 (22, 3) tile bottom edge
83.75|53.36|Barnaul
90.00|60.00|a point on L3 (22, 3) tile top edge
world-ll $ v.in.ascii input=points-ll.dat \
               columns='x double precision, y double precision, name varchar(40)' \
               output=marked_points
...
world-ll $ v.out.ascii input=marked_points
80|50|1
83.75|53.36|2
90|60|3
world-ll $

  The points on the edges of an L3 tile are to allow the `v.proj'
  result to be easily checked afterwards. (Rows of the L3 tiles
  are the stripes between consequent 10 degrees parallels, so row
  #3 is bounded by 60dN and 50dN.) The ``Barnaul'' point is to be
  checked visually.

  Then, I reproject the map:

modis-sinu $ v.proj input=marked_points location=world-ll mapset=ivan
Input Projection Parameters: +proj=latlong +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000
Input Unit Factor: 1

Output Projection Parameters: +proj=sinu +R=6371007.181 +lat_0=0.0000000000 +lon_0=0.0000000000 +no_defs
Output Unit Factor: 1
...
modis-sinu $ v.out.ascii input=marked_points
5740503.95315973|5538668.85235319|1
5581642.24347567|5912856.23809764|2
5029005.65975492|6653142.01247702|3
modis-sinu $

  The result seems not to be correct. The L3 tiles in sinusoidal
  projection are the squares comprised of 1200 x 1200 pixels each,
  and each pixel is a 926.62543306 meters square. There are 36 x
  18 L3 tiles covering the whole Earth surface.

  Rows of the tiles are numbered north-to-south, with 90..80dN row
  numbered 0th. So, 3rd row is 60..50dN, and it corresponds to
  Y-range of:

(- (* 926.62543306 1200 9) (* 926.62543306 1200 3))
;; 60dN
=> 6671703.118031999
(- (* 926.62543306 1200 9) (* 926.62543306 1200 4))
;; 50dN
=> 5559752.598359999

  While the `v.proj' results in the points being shifted by the
  amount of 19..21 km:

(- 6653142.01247702 6671703.118031999)
;; 60dN
=> -18561.105554979295
(- 5538668.85235319 5559752.598359999)
;; 50dN
=> -21083.746006809175

  When I use `proj' to ``reproject'' the data, it seems to work
  just fine:

modis-sinu $ sed -e 's,^\([^|]*\)|\([^|]*\)\(|.*\),\1 \2 \3,' \
                 < points-ll.dat \
                 | proj $(g.proj -jf) \
                 | sed -e 's,^\([^[:blank:]]*\)[[:blank:]]\+\([^[:blank:]]*\)[[:blank:]]\+\(|.*\),\1|\2\3,' \
                 | tee points-sinu.dat
5717984.13|5559752.60|a point on L3 (22, 3) tile bottom edge
5557613.28|5933367.97|Barnaul
5003777.34|6671703.12|a point on L3 (22, 3) tile top edge
modis-sinu $ v.in.ascii input=points-sinu.dat \
                 columns='x double precision, y double precision, name varchar(40)' \
                 output=marked_points_proj
...
modis-sinu $

  The points are visually at the top and at the bottom of the
  imported raster, and ``Barnaul'' seems to be at the right place.