Tutorial 3: Gridding 500 m Snow Cover Data over Colorado Using mod10_l2.pl
Requirements
Suppose we want to put some MODIS 500 m snow cover swath data covering
Colorado into an Azimuthal Equal Area spherical projection centered on
40 N 105 W with the vertical axis of the grid through the center pointing
due north towards the top of the grid. We want the grid resolution to be
500 m and we want to use a spherical earth radius value of 6370.997 km.
We want the upper left corner of the grid to be at about 44 N 110 W and
the lower right corner of the grid to be at about 36 N 100 W. We want to
grid all available MOD10_L2
"channels." These include channel 1 Snow Cover and channel 2 Snow Cover
Pixel QA. We could simply order MOD10_L2
data; however, as we can see in Supported
Data Sets, the latlon data are stored at only 5 km resolution in MOD10_L2
granules. We could have mod10_l2.pl work
with only the MOD10_L2
granules (i.e. by setting latlonlistfile to "none"); however we decide
to order the corresponding MOD03
granules as well in order to minimize interpolation error.
NOTE: To run this example, you'll need a machine with at least 100 MB
of memory and at least 300 MB of free disk space.
Searching for the Data
Let's assume that we happen to know that the afternoon of September 25,
2000 was fairly clear over Colorado and that there was a significant amount
of snow on the ground, so we use the EDG
to order two MOD10_L2
granules and the corresponding two MOD03
granules acquired on September 25, 2000 at 1745 and 1750 that appear to
cover Colorado and the neighboring region. You can use the following values
for performing the search using the EDG:
Data Set
MODIS/TERRA SNOW COVER 5-MIN L2 SWATH 500M
MODIS/TERRA GEOLOCATION FIELDS 5-MIN L1A SWATH 1KM
Search Area
Type in Lat/Lon Range:
Northern latitude: 44.0000
Southern latitude: 36.0000
Western longitude: -110.0000
Eastern longitude: -100.0000
Start Date: 2000-09-25 Time (UTC): 17:00:00
End Date: 2000-09-25 Time (UTC): 18:00:00
The search should find four granules having the following names:
MOD03.A2000269.1745.002.2000282120258.hdf
MOD03.A2000269.1750.002.2000282120547.hdf
MOD10_L2.A2000269.1745.002.2000290231013.hdf
MOD10_L2.A2000269.1750.002.2000290231230.hdf
Note that September 25, 2000 is day-of-year 269.
Ordering and Downloading the Data
Order and download the above files to some directory we'll call the tutorial_3
directory where you have at least 300 MB of free disk space. Note that
you can also download the *.met files that accompany the *.hdf files, but
the MS2GT software doesn't use them.
Creating the mod10_l2.pl Command File
Create a text file in the tutorial_3 directory called colo_2000269_1745.csh
containing the following line:
mod10_l2.pl . colo_2000269_1745 listfile.txt Colorado500.gpd 12
latlonlist.txt
This command specifies the following information (see mod10_l2.pl):
-
dirinout is "." indicating that the current directory in effect when colo_2000269_1745.csh
is invoked will contain the input and output files.
-
tag is "colo_2000269_1745" indicating that all output filenames containing
gridded data created by mod10_l2.pl will
begin with the string "colo_2000269_1745".
-
listfile is "listfile.txt" containing a list of the MOD10_L2
files to be processed (see Creating the listfile).
-
gpdfile is "Colorado500.gpd" containing a specification of the grid and
its associated map projection to use in gridding the data (see Creating
the gpd and mpp files).
-
chanlist is "12" specifying that both channels 1 (Snow Cover) and 2 (Snow
Cover PixelQA) should be gridded.
-
latlonlistfile is "latlonlist.txt" containing a list of the MOD03
files whose latitude and longitude data should be used in place of the
latlon data in the corresponding MOD10_L2
files in listfile.txt (See Creating the latlonlistfile).
-
keep is not specified, so the default value of "0" is used indicating that
intermediate chan, lat, lon, col, and row files should be deleted.
-
rind is not specified, so the default value of "50" is used. If you see
holes in the final grid that seem to correspond to the boundaries between
adjacent swath granules, then you might try increasing the rind value.
Make colo_2000269_1745.csh executable by typing:
chmod +x colo_2000269_1745.csh
Creating the listfile
Create a text file called listfile.txt in the tutorial_3 directory containing
the following two lines:
MOD10_L2.A2000269.1745.002.2000290231013.hdf
MOD10_L2.A2000269.1750.002.2000290231230.hdf
Note that we list the MOD10_L2
files to be gridded.
Creating the gpd and mpp files
See Points,
Pixels, Grids, and Cells for a description of the gpd and mpp file
formats used by the mapx library in defining a grid and its associated
map projection. In the previous two tutorials, we modified existing gpd
files and used an existing mpp file (namely N200correct.mpp)
that specified EASE-Grid projections. Here we are going to create gpd and
mpp files from scratch. We'll start with creating the mpp file, which we'll
call Colorado.mpp, in the ms2gt/grids directory (or, if you don't want
to type it in, copy Colorado.mpp from the ms2gt/tutorial_3 directory to
the ms2gt/grids directory):
Azimuthal Equal-Area
40.0 -105.0 lat0 lon0
0.0
rotation
100.0
scale (km/map unit)
40.0 -105.0 center lat lon
20.0 60.0 lat min max
-135.0 -75.00 lon min max
10.00 15.00 grid
0.00 0.00 label lat lon
1 0 0
cil bdy riv
6370.997 Earth equatorial
radius (km) -- gctp
-
The first line specifies the projection we wish to use, namely Azimuthal
Equal-Area. Since we didn't specify ellipsoid, a spherical projection is
assumed.
-
The second line specifies the map projection origin, namely 40 N 105 W.
-
The third line specifies the rotation, namely 0 degrees. This will produce
a map with the vertical axis through the center pointing due north towards
the top of the map.
-
The fourth line specifies an arbitrary scale for the map as opposed to
the grid, which will be defined by the gpd file as grid cells per map unit.
Here we define a map unit to be 100 km.
-
The fifth line specifies the center of the map which is usually (but not
necessarily) the map projection origin. Here we simply set it equal to
the map projection origin, namely 40 N 105 W.
-
The next five lines (the sixth through tenth lines) specify parameters
that would be useful to programs that produce graphic overlays (see Points,
Pixels, Grids, and Cells). They are not used by the MS2GT software,
but they need to be present in the mpp file as place holders.
-
The eleventh line, if present, specifies the equatorial radius to use instead
of the default 6371.228 km. Here we specify 6370.997 km, which happens
to be the default radius used by the General
Cartographic Transformation Package (GCTP).
In preparing the gpd file which will define our grid, we will need to know
following:
-
The name of the mpp file which will define our map, namely Colorado.mpp.
-
The number of columns and rows. We don't know these yet, but we do know:
-
There are about 40000 km / 360 degrees = 111 km/deg in longitude at the
equator or in latitude anywhere.
-
We want our grid to have 500 m per cell = 0.5 km/cell.
-
We need to span about 10 degrees in longitude (110 W to 100 W) at about
40 N. This works out to about 10 deg * cos(40 deg) * 111 km/deg / (0.5
km/cell) = 1700 cells in longitude = 1700 columns.
-
We need to span about 8 degrees in latitude (44 N to 36 N). This works
out to about 8 deg * 111 km/deg / (0.5 km/cell) = 1776 cells in latitude
= 1776 rows.
-
The number of grid cells per map unit. This is equal to 100 km/map unit
/ (0.5 km/cell) = 200 cells/map unit
-
The grid cell coordinates of the center of the map. Since we want the center
of the map to be the center of our grid, and since we have an even number
of rows and columns in our grid, the center of our grid will be at the
joint corner of the center four cells. But since grid cell coordinates
are measured relative to the center of the upper left cell, it turns out
that the grid cell coordinates of the center of the map are (1700 columns
/ 2) - 0.5 = column 849.5 and (1776 rows / 2) - 0.5 = row 887.5.
We now have all the information we need to create the gpd file which we'll
call Colorado500.gpd in the ms2gt/grids directory (or, if you don't want
to type it in, copy Colorado500.gpd from the ms2gt/tutorial_3 directory
to the ms2gt/grids directory):
Colorado.mpp map projection parameters
# azimuthal equal area
1700 1776 columns rows
#
200
grid cells per map unit
# 500 meters
849.5 887.5 origin column, row
# center at 40N 150W
Once Colorado500.gpd has been created in the ms2gt/grids directory,
we can use gtest to check that the latitude and longitude values of the
upper left and lower right corners are close to what we want and that the
center is exactly what we want:
gtest
enter .gpd file name: Colorado500.gpd
> assuming old style fixed format file
gpd: Colorado500.gpd
mpp:Colorado.mpp
forward_grid:
enter lat lon:
inverse_grid:
enter r s: 0 0
lat,lon = 43.875530 -110.300537 status = 1
col,row = -0.000061 0.000366 status = 1
enter r s: 1699 1775
lat,lon = 35.909435 -100.283951 status = 1
col,row = 1699.000000 1774.999634 status = 1
enter r s: 849.5 887.5
lat,lon = 40.000000 -105.000000 status = 1
col,row = 849.500000 887.500000 status = 1
enter r s:
enter .gpd file name:
So we see that the upper left corner values of 43.875530 -110.300537
are close to our target values of 44 N and 110 W and that the lower right
corner values of 35.909435 -100.283951 are close to our target values of
36 N and 100 W. The center values are exactly equal to our target value
of 40 N 105 W.
Creating the latlonlistfile
Create a text file in the tutorial_3 directory called latlonlist.txt containing
the following two lines:
MOD03.A2000269.1745.002.2000282120258.hdf
MOD03.A2000269.1750.002.2000282120547.hdf
Note that we list the MOD03
files from which we wish to extract the 1 km latitude and longitude data
to be used in place of the latlon data in the corresponding MOD10_L2
files in listfile.
Running the mod10_l2.pl Command File
Run the shell script containing the mod10_l2.pl
command by changing to the tutorial_3 directory, and then typing:
colo_2000269_1745.csh
You'll see lots of messages displayed while the mod10_l2.pl
script runs various IDL and C programs. In this example, the programs are:
-
extract_latlon - an IDL procedure for extracting latitude and longitude
data from a MOD10_L2
or MOD03
file. This program calls another IDL procedure, modis_ancillary_read. In
this example, extract_latlon is called twice, once for each of the two
MOD03
files. Two binary floating-point files are created per call containing
latitude and longitude data, respectively. The mod10_l2.pl
script concatenates the two latitude files and the two longitude files
to create a single latitude file and a single longitude file, and the pre-concatenated
files are deleted.
-
extract_chan - an IDL procedure for extracting channel data from a MOD10_L2
file. This program calls another IDL procedure, modis_snow_read. In this
example, extract_chan is called four times: two times for each of the two
MOD10_L2
files; on each call, channel 1 or channel 2 is extracted. One binary byte
file is created per call containing the channel data. The mod10_l2.pl
script concatenates each pair of channel files, creates two concatenated
channel files, and then deletes the pre-concatenated channel files.
-
ll2cr - a C program for converting latitude, longitude pairs to column,
row pairs for a particular grid. The grid specified in this example is
Colorado500.gpd. The concatenated latitude and longitude files are read
and two binary floating-point files are created containing column and row
numbers, respectively. The mod10_l2.pl
script then deletes the concatenated latitude and longitude files.
-
interp_colrow - an IDL procedure for interpolating column, row pairs from
a lower resolution swath format to a higher resolution swath format, in
this case from 1 km to 500 m. The interpolation must be performed on a
scan's worth of data at a time because the column and row numbers have
discontinuities at scan boundaries. The interp_colrow procedure calls a
function called congridx for each scan's worth of column and row arrays.
The congridx function is called once for the column array and once for
the row array. The congridx function first performs an extrapolation of
the given array to a slightly expanded array, which it then interpolates
(bicubic interpolation is used here) to a fully expanded array. The final
array is extracted from the fully expanded array. The mod10_l2.pl
script then deletes the pre-interpolated column and row files.
-
fornav - a C program for performing forward navigation from a swath to
a grid. In this example, fornav is called two times, once for each of the
two concatenated channel files. On each call, the column and row files
are read as well. An elliptical weighted maximum algorithm is applied during
forward navigation to minimize holes and aliasing in the gridded data.
One binary byte file is created per call containing the gridded data. The
mod10_l2.pl
script then deletes the concatenated channel files as well as the column
and row files.
The final message should contain the string:
MOD10_L2: MESSAGE: done
Examining the Results
Enter the command:
ls -l *.img
You should see something like this:
-rw-r--r-- 1 haran nsidc
3019200 Apr 27 13:04 colo_2000269_1745_rawm_snow_01700_01776.img
-rw-r--r-- 1 haran nsidc
3019200 Apr 27 13:04 colo_2000269_1745_rawm_snqa_01700_01776.img
Each file contains a gridded array of 1700 columns and 1776 rows of
binary byte values (1700 * 1776 * 1 = 3019200 bytes).
The file naming convention for gridded MOD10_L2
files is as follows:
<tag>_<conversion><weight_type>_<chan>_<columns>_<rows>.img
-
<tag> is the mod10_l2.pl tag parameter
-
<conversion> is:
-
raw - raw (1-byte unsigned integers)
-
<weight_type> is:
-
<chan> is the channel name and is one of:
-
snow - channel 1 - Snow Cover - 8-bit unsigned
-
snqa - channel 2 - Snow Cover PixelQA - 8-bit unsigned
-
<columns> is the number of columns in the grid
-
<rows> is the number of rows in the grid
Last updated: May 31, 2001 by
Terry Haran
NSIDC-CIRES
449 UCB
University of Colorado
Boulder, CO 80309-0449
303-492-1847
tharan@colorado.edu