Working with MODIS L1B from scratch. 2. subsetting and reprojecting

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 that are shown here can also be made with the "gdalinfo" and "gdalwarp" tools that are part of the GDAL suite. (I will add an update on this later!) - We begin with the two hdf files downloaded in the previous step: MOD021KM.A2013187.1035.005.2013187194950.hdf MOD03.A2013187.1035.005.2013187174113.hdf 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 are describes as: OBJECT=DataField_5 DataFieldName="EV_250_Aggr1km_RefSB" DataType=DFNT_UINT16 DimList=("Band_250M","10*nscans","Max_EV_frames") END_OBJECT=DataField_5 OBJECT=DataField_8 DataFieldName="EV_500_Aggr1km_RefSB" DataType=DFNT_UINT16 DimList=("Band_500M","10*nscans","Max_EV_frames") END_OBJECT=DataField_8 - 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: 65.00N 8.00E 55.00N 22.00E 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 65.00N 8.00E 55.00N 22.00E ) - 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: swath2grid -pf=parameter.prm 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 ------------------------------------------------------------------ debug 1 debug 2 General processing info ----------------------- input_filename: MOD021KM.A2013187.1035.005.2013187194950.hdf geoloc_filename: MOD03.A2013187.1035.005.2013187174113.hdf output_filename: sweden output_filetype: GEOTIFF output_projection_type: Universal Transverse Mercator (UTM) output_zone_code: 33 output_ellipsoid: GRS 1980/WGS 84 output_datum: WGS84 resampling_type: NN 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: sweden_EV_250_Aggr1km_RefSB_b0.tif sweden_EV_250_Aggr1km_RefSB_b1.tif sweden_EV_500_Aggr1km_RefSB_b0.tif sweden_EV_500_Aggr1km_RefSB_b1.tif sweden_EV_500_Aggr1km_RefSB_b2.tif sweden_EV_500_Aggr1km_RefSB_b3.tif sweden_EV_500_Aggr1km_RefSB_b4.tif 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.