Col-Row Subsetting of MOD29P1D Granules V0.2

Contents

Version History
Purpose
MOD29P1D Tiles
Converting Local Grid Coordinates to Absolute Grid Coordinates
The X-Y Coordinate System
Converting Absolute Grid Coordinates to X-Y Coordinates
Converting X-Y Coordinates to Local Grid Coordinates
Example

Version History

V0.2 released February 25, 2004.

V0.1 released February 12, 2004.

V0.0 released July 9, 2002. 

Purpose

The purpose of this document is to describe a plan describing how col-row subsetting between EDG/ECS and HEW might work using MOD29P1D granules (tiles) as an example.

The intended audience for this document is anyone interested in defining how col-row subsetting in the EDG/ECS/HEW environment might work.

MOD29P1D Tiles

See MODLAND Polar Grids for a description of the MODLAND Polar Grids used for MOD29P1D MODIS/TERRA SEA ICE EXTENT DAILY L3 GLOBAL 1KM EASE-GRID DAY tiles.

Each hemisphere consists of an array of 19x19 tiles numbered 0 to 18, with the "h" tile numbers listed across the top and bottom, and the "v" tile numbers listed down the sides. For the northern hemisphere, the origin of the projection is at latitude 90, longitude 0, and the v tile numbers are in the range 0 to 18. For the southern hemisphere, the origin is at latitude -90, longitude 0, and v tile numbers are in the range 20 to 38.

For the MOD29P1D product, each tile is 951x951, so the total grid for each hemisphere is 18069x18069. Let's adopt a col-row numbering system with the origin at the center of the upper left pixel. Column values increase to the right, and row values increase downward. The column and row values of the center of each pixel within a tile (the local grid) are in the range 0 to 950.  The column and row values of the center of each pixel in the absolute grid are in the range 0 to 18068, with the center of the projection in the center of the pixel having absolute (col, row) coordinates of (9034, 9034).
 

Converting Local Grid Coordinates to Absolute Grid Coordinates

We can convert local grid coordinates to absolute grid coordinates using:

abs_col = h * cols_per_tile + loc_col
abs_row = (v - v_offset) * rows_per_tile + loc_row

where
  h is the "h" tile number
  v is the "v" tile number
  loc_col is the local column number
  loc_row is the local row number
  v_offset = 0 for the northern hemisphere
                 = 20 for the southern hemisphere
  cols_per_tile = 951 (XDim in the structural metadata in a MOD29P1D tile)
  rows_per_tile = 951 (YDim in the structural metadata in a MOD29P1D tile)

See loc2abs.pro for an IDL program that converts local grid coordinates to absolute grid coordinates using the above algorithm.

Converting Absolute Grid Coordinates to Local Grid Coordinates

We can convert absolute grid coordinates to local grid coordinates using:

h = fix(abs_col / cols_per_tile)
v = fix(abs_row / rows_per_tile) + v_offset
loc_col = abs_col mod cols_per_tile
loc_row = abs_row mod rows_per_tile

See abs2loc.pro for an IDL program that converts absolute grid coordinates to local grid coordinates using the above algorithm.

The X-Y Coordinate System

The x-y coordinate system, which is used to denote the upper left and lower right corners within an HDF-EOS file, is in units of meters, with the origin at the center of the projection, x increasing to the right, and y increasing upward. When denoting upper left and lower right x-y coordinates, the convention used by the HDF-EOS library is to use the upper left corner of the upper left pixel and the lower right corner of the lower right pixel, respectively. The x and y resolution for a tile can be calculated using:

meters_per_pixel_x = (tile_lr_x - tile_ul_x) / cols_per_tile
meters_per_pixel_y = (tile_ul_y - tile_lr_y) / rows_per_tile

where
  tile_ul is the two element array UpperLeftPointMtrs(x,y) in the structural
     metadata in a MOD29P1D tile.
  tile_lr is the two element array LowerRightMtrs(x,y) in the structural
     metadata in a MOD29P1D tile.

The meters_per_pixel_x and meters_per_pixel_y values should be the same for all tiles in a dataset.

Converting Absolute Grid Coordinates to X-Y Coordinates

Let's assume that EDG/ECS is going to obtain upper left and lower right absolute grid col-row pairs from the user defining a subset region for a number of MOD29P1D tiles for a single known hemisphere (EDG/ECS can deduce the hemisphere from the v values in the MOD29P1D filenames) and will define this absolute col-row subset region using the following variable names:

subset_ul_abs_col
subset_ul_abs_row
subset_lr_abs_col
subset_lr_abs_row

Let's further assume that EDG/ECS wants to convert the absolute col-row pairs to x-y pairs for the purpose of sending the x-y pairs to HEW in order to define the subset region. Then EDG/ECS needs to know that for MOD29P1D, the center of the projection is:

center_abs_col = 9034.0
center_abs_row = 9034.0

and the resolution of the grid in meters per pixel is:

meters_per_pixel_x = 1002.701
meters_per_pixel_y = 1002.701

Then EDG/ECS can compute the upper left and lower right x-y pairs for the upper left and lower right corners of the upper left and lower right pixels, respectively, in the subset region:

subset_ul_x = (subset_ul_abs_col - center_abs_col - 0.5) * meters_per_pixel_x
subset_ul_y = (center_abs_row - subset_ul_abs_row + 0.5) * meters_per_pixel_y
subset_lr_x = (subset_lr_abs_col - center_abs_col + 0.5) * meters_per_pixel_x
subset_lr_y = (center_abs_row - subset_lr_abs_row - 0.5) * meters_per_pixel_y

See abs_colrow2xy.pro for an IDL program that converts an absolute col-row subset region to an x-y subset region using the above algorithm.

Converting X-Y Coordinates to Local Grid Coordinates

Once HEW gets the x-y pairs defining the subset region for all the tiles, it needs to compute the upper left and lower right local col and row values defining the subset region for each individual tile. From the structural metadata in each tile, HEW knows (see above):

(tile_ul_x, tile_ul_y) = UpperLeftPointMtrs(x,y)
(tile_lr_x, tile_lr_y) = LowerRightMtrs(x,y)
cols_per_tile = XDim
rows_per_tile = YDim

Then for each tile, HEW computes:

meters_per_pixel_x = (tile_lr_x - tile_ul_x) / cols_per_tile
meters_per_pixel_y = (tile_ul_y - tile_lr_y) / rows_per_tile
subset_ul_loc_col = max(round((subset_ul_x - tile_ul_x) / meters_per_pixel_x), 0)
subset_ul_loc_row = max(round((tile_ul_y - subset_ul_y) / meters_per_pixel_y), 0)
subset_lr_loc_col = min(round((subset_lr_x - tile_ul_x) / meters_per_pixel_x) - 1, cols_per_tile - 1)
subset_lr_loc_row = min(round((tile_ul_y - subset_lr_y) / meters_per_pixel_y) - 1, rows_per_tile - 1)
if ((subset_ul_loc_col < cols_per_tile) and
     (subset_ul_loc_row < rows_per_tile) and
     (subset_lr_loc_col >= 0) and
     (subset_lr_loc_row >= 0)) then
        subset_this_tile = TRUE
else
        subset_this_tile = FALSE

See xy2loc_colrow.pro for an IDL program that converts an x-y subset region to a local col-row subset region for a particular tile using the above algorithm.

Example

Suppose we want to order and subset some MOD29P1D data for the Beaufort Sea for April 1, 2002. We perform the following steps (steps 3, 4, and 6 are performed by the IDL script example.pro):
Terry Haran <tharan@nsidc.org>

Last modified: Thu Feb  25 14:32:51 MST 2004