In the second part of this series, we will see how to extract subsets of the hdf files we downloaded and reproject them. We will use the programs “dumpmeta” and “swath2grid”, that are part of the MODIS Reprojection Tool Swath suite. These can be downloaded from the LP DAAC site. The package is free of cost, but you will have to register to gain access to it. Instructions for installing it come with the package. [In Linux this is is just a matter of copying the uncompressed tree into a suitable directory, like /usr/local/MRTSwath, and then adding the ‘bin’ subfolder to the $PATH]. All procedures can also be performed with the “gdalinfo” and “gdalwarp” tools that are part of the GDAL suite.
– We begin with the two hdf files downloaded in the previous step:
The first one is the image data, while the second includes the ancillary information on georeference, sun position, etc.
– We then start by exploring the contents of our file using the program “dumpmeta”:
dumpmeta MOD021KM.A2013187.1035.005.2013187194950.hdf metadata.dat
We open the output file in a text editor and look for the naming of the data granules. Under the heading “GROUP=DataField” we will find the names. In this example we look for bands acquired in the visible and mid-infrared regions of the spectrum: bands 1-7. These bands are originally acquired using 250 m (bands 1 and 2, red and near-infrared) and 500 m pixel size (bands 3-7, from the blue to the mid-infrared). In the 1KM product we downloaded, these bands have been resampled. The names look like this: “EV_250_Aggr1km_RefSB” and “EV_500_Aggr1km_RefSB”, that is, 250m and 500m that have been aggregated into 1km bands.
In the metadata file this is described as:
– We then select the coordinates of the upper left and lower right corner we want to cut from the image. In our case, we will use:
but we want them reprojected in a suitable UTM projection, in our case UTM33, so using the coordinates are approximately:
65.00N 8.00E 170432.16 7226731.39 0.00
55.00N 22.00E 947397.04 6117225.18 0.00
This was done using the program “cs2cs” that comes with the package proj4:
cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs -r -E
– Now we have all elements we need to extract a subset. We put all relevant information in a parameter file “parameter.prm”. The contents of this file are as follows:
INPUT_FILENAME = MOD021KM.A2013187.1035.005.2013187194950.hdf
OUTPUT_FILENAME = sweden
OUTPUT_FILE_FORMAT = GEOTIFF_FMT
GEOLOCATION_FILENAME = MOD03.A2013187.1035.005.2013187174113.hdf
INPUT_SDS_NAME = EV_250_Aggr1km_RefSB; EV_500_Aggr1km_RefSB
KERNEL_TYPE = NN
OUTPUT_PROJECTION_NUMBER = UTM
OUTPUT_PROJECTION_PARAMETER = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
OUTPUT_SPACE_UPPER_LEFT_CORNER = 170000 7230000
OUTPUT_SPACE_LOWER_RIGHT_CORNER = 950000 6120000
OUTPUT_SPATIAL_SUBSET_TYPE = PROJ_COORDS
OUTPUT_PROJECTION_SPHERE = 8
OUTPUT_PROJECTION_ZONE = 33
OUTPUT_PIXEL_SIZE = 1000
In order, these lines specify the input image filename, some generic output name, the output file format, the geolocation file name, the data granules to be exported, the resampling method (here nearest neighbours), the cartographic projection, a list of parameters related to the projection, the coordinates of the UL corner, same with the LR corner, a label specifying that we have projected coordinates, a projection sphere (8 = GRS80/WGS84), an UTM zone and the output pixel size in meters.
– We then execute the command:
which dumps the following output in the screen:
MODIS Reprojection Tool Swath (v2.2 September 2010)
Start Time: Fri Nov 15 17:28:28 2013
General processing info
output_projection_type: Universal Transverse Mercator (UTM)
output_ellipsoid: GRS 1980/WGS 84
output projection parameters: 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
SDS name #bands in SDS #bands to process
1) EV_250_Aggr1km_RefSB 2 2
2) EV_500_Aggr1km_RefSB 5 5
Processing EV_250_Aggr1km_RefSB, 0 ...
output lines/samples: 1110 780
output pixel size: 1000.0000
output data type: same as input
output upper left corner: lat 65.02868042 long 7.98322664
output lower right corner: lat 55.03222632 long 22.03079051
input lines/samples: 2030 1354
input resolution: 1 km
%% complete: 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
… and so on until is finished with all 7 bands.
– We finally have 7 new image files in our folder:
These files contain the first 7 bands of the reprojected subset of the image data we downloaded. The image type is 16 bit unsigned integer. The naming b0, b1, etc is related to the data granules. In the 250_Aggr1km case, b0 and b1 mean bands 1 (red) and 2 (near infrared). The band specifications can be found here.
We can combine the images sweden_EV_250_Aggr1km_RefSB_b0.tif, sweden_EV_500_Aggr1km_RefSB_b1.tif and sweden_EV_500_Aggr1km_RefSB_b0.tif in RGB channels to obtain a true color composite. We do this using the program “gdal_merge.py” that comes with GDAL:
gdal_merge.py -separate -o sweden_2013187_143.tif sweden_EV_250_Aggr1km_RefSB_b0.tif sweden_EV_500_Aggr1km_RefSB_b1.tif sweden_EV_500_Aggr1km_RefSB_b0.tif
The results look like this:
The image in the example above has been actually reduced and enhanced to better fit into this page. This was done simply using ImageMagick:
convert sweden_EV_250_Aggr1km_RefSB_b0.tif sweden_EV_500_Aggr1km_RefSB_b1.tif sweden_EV_500_Aggr1km_RefSB_b0.tif -channel RGB -combine -geometry 25% -alpha off -equalize sweden_2013187_143.png
Of course, all these procedures are very tedious to do by hand and should be automatized for large scale projects in which hundreds of images are involved. In such a case, all of the above must be embedded in a script, using for example Perl or Python.
In the next post, we will see how to convert these images to radiance at the top of the atmosphere.
Note: originally published in my old blog Tales of Ice and Stone on 2013-11-15