[R-sig-Geo] Pixel reference location?

[Lyndon Estes]

Dear R-Sig-Geo,

I have a question regarding pixel referencing with the raster package,
from which I am using the crop function to resize MODIS imagery.  I
was under the impression that raster refers to the upper left of the
upper left pixel in defining the spatial extent, but I am getting
results that seem to suggest that the center of the pixel is being
used.

What I am doing: I have written a function that calls the MODIS
HEGtools to read-in, stitch, and reproject the image into an Albers
Equal Area Conic projection.  The HEGtool has the capacity to subset
the image, but the resulting subset image is not well rectified with
the result one gets when stitching two full size images together
(probably because of how I have defined projection parameters), so I
am instead adding an extra step that resizes the image with the crop
function, using a subset of the image that I created using ENVI.

Here are the values for the crop image:

crop.img
class       : RasterLayer
filename    : MOD250m_extenttemplate
nrow        : 6116
ncol        : 4200
ncell       : 25687200
min value   : ?
max value   : ?
projection  : +proj=aea +lat_1=-18 +lat_2=-32 +lat_0=-30 +lon_0=24
+x_0=0 +y_0=0 +ellps=clrk66 +datum=NAD27 +units=m +no_defs
+nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat
xmin        : -120863.9
xmax        : 852092.8
ymin        : -497022.6
ymax        : 919787.6
xres        : 231.6564
yres        : 231.6564

(Note the seemingly bizarre projection parameters, which is another
story, and is created by the HEGtool, which is supposed to use WGS4,
but instead seems to spit out NAD27).

Now, after the crop function is run, my resulting cropped MODIS image
ends up with the following parameters:

raster(out.name)
class       : RasterLayer
filename    : MOD13Q1.A2005225.h20v11.mosaic4.NDVI.tif
nrow        : 6116
ncol        : 4200
ncell       : 25687200
min value   : ?
max value   : ?
projection  : +proj=aea +lat_1=-18 +lat_2=-32 +lat_0=-30 +lon_0=24
+x_0=0 +y_0=0 +ellps=clrk66 +datum=NAD27 +units=m +no_defs
+nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat
xmin        : -120748.1
xmax        : 852208.6
ymin        : -496906.8
ymax        : 919903.5
xres        : 231.6564
yres        : 231.6564

Looking at the xmins of the two images, -120863.9 - -120748.1 =
-115.8, which is equal to half a pixel.  Same for the ymaxs.

So, is this an issue of where in the pixel raster assigns the UL
reference coordinate, and, if so, is there a way to specify how the
reference is being done?  If not, am I doing something else wrong?

Thanks in advance for any advice (function code follows).

Cheers, Lyndon


ModStitch <- function(modpath, hegpath, modnames, cropimage) {
# Uses the HEGtool from MODIS LPDAAC to read in, stitch, and reproject
the MODIS NDVI and VI quality bands
# References: http://spatial-analyst.net/book/MODIS_download_HR.R
# Args:
#   modpath: Path to directory holding MODIS data
#   hegpath: Path to the HEGtool installation
#   modnames: Vector containing MODIS image names to process
#   crop image: A raster defining the extent to which the images
should be resized
# Returns:
#   Separate stitched MODIS NDVI and VI Quality bands, reprojected to
Albers Equal Area conic (for South
#   Africa), saved as a GEOtiff

  # Set tool and environment variables
  hegtool         <- paste(hegpath, "/heg/bin/hegtool", sep = "")
  grd.stitch.tool <- paste(hegpath, "/heg/bin/subset_stitch_grid", sep = "")
  pgs.loc         <- paste(hegpath, "/heg/TOOLKIT_MTD", sep = "")
  mrt.data.loc    <- paste(hegpath, "/heg/data", sep = "")
  Sys.setenv(PGSHOME = pgs.loc)  # Environment variable for HEGtool
  Sys.setenv(MRTDATADIR = mrt.data.loc)  # Environment variable for HEGtool

  # Set path to MODIS data directory
  setwd(modpath)

  # And process the specified MODIS images using hegtool to capture
header file information
  mod.list <- vector("list", length(modnames))  # List to capture
header file output

  # Call crop image
  crop.img <- raster(cropimage)

  # Create and capture HegHdr.hdr files
  for (i in 1:length(mod.list)) {
      system(paste(hegtool, "-h", modnames[i]), intern = FALSE, wait = TRUE)
      mod.list[[i]] <- scan("HegHdr.hdr", what = "character", sep = "\n",)
      names(mod.list)[i] <- modnames[i]
  }

  # Then stitch and reproject both the NDVI layers and the error flag layers

  # Set up a couple of variables
  which.var <- c("250m 16 days NDVI|", "250m 16 days VI Quality|")
  which.abb <- c("NDVI", "qual")

  for(k in 1:length(which.var)) {  # Loop through NDVI, then VI

    # First write the stitch parameters file

    # Stitch parameter file name
    pr.name  <- paste("stitch_", k, ".prm", sep = "")  # 1 = stitch
for NDVI, 2 for VI quality
    out.name <- paste(substr(modnames[1], 1, 23), ".mosaic4.",
which.abb[k],".tif", sep = "")

    # Open connection to write stitch file
    filename <- file(pr.name, open = "wt")
    write("", filename, append = TRUE)
    write("NUM_RUNS = 1", filename, append = TRUE)  # Not sure if this
needs to be a variable
    write("", filename, append = TRUE)
    write("BEGIN", filename, append = TRUE)
    write("NUMBER_INPUTFILES = 2",filename, append = TRUE)  # Alter
this if more than 1 file is to be stitched
    write(paste("INPUT_FILENAMES = ", modnames[1], "|", modnames[2],
sep = ""), filename, append = TRUE)
    write("OBJECT_NAME = MODIS_Grid_16DAY_250m_500m_VI|", filename,
append = TRUE)
    write(paste("FIELD_NAME = ", which.var[k], sep = ""), filename,
append = TRUE)
    write("BAND_NUMBER = 1", filename, append = TRUE)

    # Pulls coordinates for upper left from first header file
    write(paste("SPATIAL_SUBSET_UL_CORNER = ", "( ",
                substr(mod.list[[1]][grep("GRID_UL_CORNER_LATLON=",
mod.list[[1]])], 23, 48), " )", sep = ""),
          filename, append = TRUE)

    # Pulls coordinates for lower right from last header file
    write(paste("SPATIAL_SUBSET_LR_CORNER = ", "( ",

substr(mod.list[[length(mod.list)]][grep("GRID_LR_CORNER_LATLON=",
                  mod.list[[length(mod.list)]])], 23, 48), " )", sep = ""),
          filename, append = TRUE)

    write("OUTPUT_OBJECT_NAME = MODIS_Grid_16DAY_250m_500m_VI|",
filename, append = TRUE)
    write(paste("OUTGRID_X_PIXELSIZE = ",
                substr(mod.list[[1]][grep("GRID_PIXEL_SIZE=",
mod.list[[1]])], 17, 26), sep = ""),
          filename, append = TRUE)
    write(paste("OUTGRID_Y_PIXELSIZE = ",
                substr(mod.list[[1]][grep("GRID_PIXEL_SIZE=",
mod.list[[1]])], 17, 26), sep = ""),
          filename, append = TRUE)
    write(paste("OUTPUT_FILENAME = ", out.name, sep = ""), filename,
append = TRUE)
    write("SAVE_STITCHED_FILE = NO", filename, append = TRUE)
    write(paste("OUTPUT_STITCHED_FILENAME = sasquatch", k, ".hdf", sep
= ""), filename, append = TRUE)
    write("OUTPUT_TYPE = GEO", filename, append = TRUE)
    write("RESAMPLING_TYPE = NN", filename, append = TRUE)
    write("OUTPUT_PROJECTION_TYPE = ALBERS", filename, append = TRUE)
    write("OUTPUT_PROJECTION_PARAMETERS = ( 0.0 0.0 -18.0 -32.0 24.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0  )",
          filename, append = TRUE)
    write("END", filename, append = TRUE)
    write("", filename, append = TRUE)
    close(filename)

    # Then call the stitch tool and process it
    stitch.com <- paste(grd.stitch.tool, "-P", pr.name)
    system(stitch.com)  # Works

    # Then use Raster function to resize because subsetting with MODIS
tool results in offset
    to.crop <- raster(out.name)
    cropped <- crop(to.crop, crop.img, filename = out.name, overwrite
= TRUE)  # Overwrites existing images
  }  # Close loop
}  # End function

ModStitch(modpath = path.ndvi.dat, hegpath = hegpath, modnames =
modnames, cropimage = envi.img)










-- 
Lyndon Estes
Research Associate
Woodrow Wilson School
Princeton University
+1-609-258-2392 (o)
+1-609-258-6082 (f)
+1-202-431-0496 (m)
lestes at princeton.edu